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PREFACE 

The  woik  reported  herein  was  performed  by  the  Arnold  Engineering  Development  Center 
(AEDC),  Air  Force  Materiel  Command  (AFMC)  under  Program  Element  921  EOS,  AEDC  Project 
Number  01 1 8,  at  the  request  of  the  AEDC  Directorate  of  Technology  (AEDCVDOT).  The  Air  Force 
Project  Manager  was  Capt.  S.  G.  Tennent.  The  results  were  obtained  by  Micro  Craft  Technology/ 
AEDC  Operations,  support  contractor  for  aerodynamic  testing  at  the  AEDC,  AFMC,  Arnold  Air 
Force  Base,  TN.  The  work  was  conducted  during  the  period  October  1,  1990  through  September 
30, 1994.  The  document  was  submitted  for  publication  on  May  5, 1995. 
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1.0  INTRODUCTION 


1.1  GENERAL 

This  work  is  motivated  by  the  need  to  understand  high-speed  continuum  flow  which  occurs 
within  the  flight  envelope  of  many  of  the  vehicles  of  current  interest,  including  single-stage-to-orbit 
vehicles,  aero-assisted  orbit-transfer  vehicles,  and  particularly  interceptor  missiles.  Personnel  at  the 
Arnold  Engineering  Development  Center  (AEDC)  are  now  engaged  in  efforts  to  build/upgiade 
various  test  facilities  to  make  them  capable  of  testing  the  new  generation  of  vehicle  systems  in  the 
high  velocity  regime  (Refs.  1-5).  These  facilities  include  the  hypersonic  tunnels,  the  ballistic 
ranges,  and  the  arc  heater  tunnels.  For  example,  the  hypersonic  tunnels  can  be  used  to  perform  jet 
interaction  and  staging  studies  for  the  interceptors.  The  ballistic  ranges  can  be  used  for  testing  the 
lethality  of  these  missiles.  The  arc  heater  tunnels  can  be  used  to  test  materials  for  interceptor  nose 
tips.  However,  it  is  important  to  emphasize  that  these  facilities  will  be  used  to  generate  conditions 
that  only  simulate  flight  conditions.  In  some  cases  there  may  be  significant  differences  between  the 
simulated  and  actual  flight  conditions.  Flight  testing  can  be  used  to  obtain  data  that  are  free  of  the 
simulation  errors  but  at  great  expense.  On  the  other  hand,  computational  simulations  are  not  subject 
to  the  limitations  specific  to  ground  testing.  Therefore,  computations  can  be  used  together  with 
ground  testing  to  accurately  capture  the  phenomena  under  investigation.  Thus  a  synergistic,  truly 
integrated  combination  of  computations,  ground  testing,  and  flight  testing  is  required  for 
hypersonic  systems  development  (Ref.  6). 

To  support  this  synergism,  an  effort  has  been  initiated  to  develop  the  computational  capability 
to  simulate  flows  in  the  test  facilities  with  high  accuracy.  Based  on  recent  scaling  studies,  it  is 
assumed  that  a  flow  solver  which  accurately  predicts  the  fluid  dynamics  of  the  test  cells  can  be  used 
to  extrapolate  to  free  flight  conditions  with  the  same  accuracy,  provided  valid  physical  models  are 
available  for  the  conditions  of  interest. 

To  establish  the  proper  context  for  the  present  effort,  note  that  the  high-speed  flows  mentioned 
above  are  characterized  by  nonequilibrium  thermochemical  processes.  Such  effects  are  signiflcant 
because  the  transition  time  for  a  fluid  particle  through  a  region  of  interest  is  shorter  than  the  time 
required  for  particle-particle  collisions  to  bring  the  gas  into  equilibrium.  Therefore,  the 
development  of  a  flow  solver  to  study  these  flows  necessarily  involves  the  modeling  of  thermo- 
chemical  nonequilibrium.  Furthermore,  the  geometric  complexities  inherent  in  the  simulation  of  the 
flows  around  the  vehicles  mentioned  require  the  use  of  domain  decomposition  techniques  which 
have  reached  maturity  in  the  chimera  methodology  (Ref.  7). 

1.2  BACKGROUND 

Computational  fluid  dynamics  has  matured  to  a  stage  where  it  is  possible  to  compute  transonic 
flow  fields  about  complex  three-dimensional  bodies  and  bodies  in  relative  motion  (Refs.  8  and  9). 
However,  there  is  a  need  to  develop  algorithms  for  complex  three-dimensional  bodies  and  bodies 
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in  relative  motion  in  hypersonic  flow  fields  which  require  reasonable  computational  resources  (Ref. 
10).  Some  important  features  of  a  good  computational  algorithm  are:  (1)  computational  efficiency 
to  minimize  the  arithmetic  operation  count  and  computer  memory  requirements,  (2)  adaptability  to 
computational  domains  consisting  of  either  blocked  or  non-structuied  grids,  (3)  fast  convergence, 
(4)  high  accuracy  which  does  not  degrade  readily  because  of  a  sensitivity  to  the  choice  of  parameter 
settings,  and  (5)  flexibility  to  allow  transport  equations  (i.e.,  species,  energy,  or  turbulence)  to  be 
added  or  deleted  with  minimal  effort. 

Solution  algorithms  can  generally  be  categorized  as  either  explicit  or  implicit  methods. 
Explicit  methods  have  low  arithmetic  operation  counts  because  they  do  not  require  matrix 
inversions  in  their  solution  procedures.  These  methods  are  limited  to  small  computational  time 
steps  for  numerical  stability.  Hence,  a  large  number  of  iterations  is  required  for  convergence  to  a 
steady-state  solution.  Implicit  methods  are  stable  for  much  larger  computational  time  steps  and 
generally  require  fewer  iterations  to  converge.  Unfortunately,  implicit  methods  require  matrix  in¬ 
versions,  which  are  computationally  intensive. 

Recently,  several  investigators  have  proposed  numerical  algorithms  which  are  globally 
explicit  and  locally  implicit.  These  algorithms  are  called  locally  implicit  methods  (LIM)  or  point- 
implicit  methods.  Locally  implicit  algorithms  have  demonstrated  the  best  features  of  both  explicit 
and  implicit  algorithms.  Reddy  and  Jacocks  (Ref.  1 1)  developed  a  two-dimensional,  finite-volume, 
locally  implicit  scheme  for  solution  of  the  Euler  equations.  The  scheme  uses  a  relaxation  method 
based  on  a  modification  to  the  one-step  Gauss-Seidel-Newton  iteration  and  does  not  require  the 
solution  of  any  matrix  equations.  The  scheme  was  applied  to  the  two-dimensional  Navier-Stokes 
equations  by  Nayani  (Ref.  1 2)  and  Towne  (Ref.  1 3).  Reddy  and  Benek  (Ref.  1 4)  developed  a  three- 
dimensional  thin-layer  Navier-Stokes  LIM  algorithm  which  incorporates  the  chimera  domain 
decomposition  procedure  developed  by  Benek  et  al.  (Ref.  7).  Hwang  and  Liu  (Ref.  15)  developed 
a  two-dimensional  finite-element  LIM  scheme.  Tramel  (Ref.  16)  developed  a  stable  shock¬ 
capturing  locally  implicit  scheme  using  a  modified  one-step  red/black  Jacobi-Newton  iteration. 
This  procedure  was  extended  to  reacting  flows  by  Tramel,  et  al.  (Ref.  17),  and  applied  to 
vibrationally  relaxing  flows  in  nozzles  by  Limbaugh  et  al.  (Ref.  18).  Bussing  and  Murman  (Ref.  19) 
used  a  point-implicit  treatment  of  the  chemical  source  terms  for  compressible  flow  problems  with 
finite-rate  chemistry.  Eberhardt  and  Imlay  (Ref.  20)  developed  a  diagonal  treatment  of  the  chemical 
source  term  Jacobian  for  the  computation  of  nonequilibrium  flow  fields.  Gnoffo  (Ref.  21) 
developed  a  finite-volume,  point-implicit  relaxation  algorithm  for  the  Navier-Stokes  equations, 
and  extended  this  algorithm  to  flows  with  chemical  and  thermal  nonequilibrium  (Ref.  22). 

1.3  OVERVIEW 

In  this  report,  a  three-dimensional,  time-accurate,  locally  implicit  algorithm  for  the  solution 
of  compressible  viscous  flow  problems  with  thermo-chemical  nonequilibrium  is  described.  The 
algorithm  is  known  as  the  nonequilibrium  diagonal  approximate  Newton’s  algorithm  (NEDANA). 
The  new  scheme  is  designed  for  vectorization  in  all  coordinate  directions.  The  scheme  incorporates 
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the  chimera  domain  decomposition  procedure.  This  allows  the  scheme  to  be  applicable  to  both 
complex  configurations  as  well  as  bodies  in  relative  motion.  The  code  also  incorporates  the  AEDC- 
developed  chemistry  package  NEQPAK  (Ref.  23).  NEQPAK  provides  the  necessary  chemical, 
thermodynamic,  and  transport  properties  that  are  required  to  simulate  the  flow  of  a  gas  in  thermo¬ 
chemical  nonequilibrium.  The  details  of  the  thermo-physical  models  and  the  numerical  method  are 
presented  in  this  report  along  with  the  results  from  a  sequence  of  test  problems. 

2.0  MODELING  EQUATIONS 
2.1  CONSERVATION  EQUATIONS 

The  NEDANA  code  solves  the  conservation  equations  which  describe  the  motion  of  a 
compressible  viscous  fluid  in  thermo-chemical  nonequilibrium  (Refs.  23,  24, 25).  A  more  general 
development  of  these  equations  is  contained  in  Appendix  A.  These  equations  are  one  mass 
conservation  equation  for  each  of  the  ns  chemical  species  present  in  the  flow. 

Op,  .  0(p,Uj)  0(p,u'l^)  _ 

dt  Oxj  ~  Ox,  (1) 


and  three  mass-averaged  momentum  equations,  (< »  1, 3), 


Ojpu,)  0(pu,Uj  +p6,j)  _  pTij 

dt  dxj  dxj '  (2) 


and  a  total  energy  equation, 

dE  d{[E  +  p)uj) 
Oi  dxj 


dxj 


0>}3 

dxj 


d 


»=i 


(3) 


where 


E  = 


El  +  -puiui, 


(4) 


and  E^is  the  internal  energy  of  the  mixture. 

The  above  equations  are  applicable  to  any  fluid  for  which  a  continuum  description  is 
appropriate.  If  the  collision  rate  among  the  individual  particles  in  the  fluid  is  high  enough  to  ensure 
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that  the  species  internal  degrees  of  freedom  are  in  thermal  equilibrium  with  the  translational  tem¬ 
perature  of  the  fluid,  then  the  above  equation  set,  plus  certain  constituent  relations  described  in  Sec. 
2.2,  provides  a  complete  description  of  the  fluid.  This  equation  set  is  referred  to  as  a  one- 
temperature  model.  However,  if  the  internal  degrees  of  freedom  are  not  in  thermal  equilibrium,  then 
additional  conservation  equations  must  be  solved.  In  the  present  work,  a  two-temperature  model  is 
assumed.  In  this  model,  the  distribution  of  the  rotational  degrees  of  freedom  is  assumed  to  be  given 
by  a  Boltzmann  distribution  whose  Boltzmann  temperature  is  equal  to  the  translational  temperature 
of  the  fluid,  while  the  vibrational/electronic  excitation  degrees  of  freedom  ate  characterized  by  a 
separate  Boltzmann  temperature  (Ref.  26).  This  leads  to  the  specification  of  an  additional 
conservation  equation  for  the  vibrational-electronic  energy,  In  the  absence  of  ionization  and 
radiation,  this  equation  takes  the  form 


d(EvUj) 
dt  dxj 


(5) 


In  this  case,  the  inlemal  energy  of  the  mixture,  ,  is  the  sum  of  a  translational/rotational  part,  E  , 
and  a  vibrational/electronic  portion,  .  i.  e., 


El  =  Elr  -f-  Ev, 


(6) 


where 


Etr  —  ^  \  PjCtr.g. 

J=1 


(7) 


and 


ns 

Ev  = 

,S=1 


(8) 


2.2  AEROTHERMAL  MODELS 

In  order  to  completely  specify  the  fluid  being  simulated,  certain  constituent  relations  are 
required:  (1)  the  thermal  equation  of  state  of  the  gas,  p  =  p(p,,A<E.7);  (2)  the  caloric  equation  of 
state  for  each  species,  h,  =  h,{T,  =  eSJ,  +  pjp^,  (3)  the  shear  stress  tensor,  t„;  (4)  the 

conductive  heat  flux,  (5)  the  diffusion  velocity  of  species  s,  uy,  (6)  the  chemical  source  term  of 
species  s,  co^;  and  (7)  the  vibrational/electronic  energy  source  term,  NEQPAK  is  used  to  provide 
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these  quantities.  NEQPAK  has  the  capability  to  apply  many  thermo-chemical  models.  Only  those 
portions  of  NEQPAK  relevant  to  the  present  work  will  be  discussed  below.  In  particular,  if  ioniza¬ 
tion  was  present  while  the  gas  was  in  a  state  of  thermal  nonequilibrium,  then  some  of  the  formulas 
presented  below  would  need  to  be  modified.  A  detailed  description  of  NEQPAK  is  found  in  Ref.  23. 

The  gas  is  assumed  to  be  composed  of  a  mixture  of  thermally  perfect  gases.  The  static  pressure 
of  the  mixture  is  then  the  sum  of  the  partial  pressures  of  the  constituent  gases  (Dalton’s  law  of 
partial  pressures), 


where 


716 

p = 

s=i 


(9) 


Pa  — 


(10) 


and  is  the  species  gas  constant.  The  species  gas  constant  is  defined  in  terms  of  the  universal  gas 
constant, 72.,  and  the  species  molecular  weight,  by  the  relationship 

R 

'  M,'  (11) 


For  a  gas  in  thermal  equilibrium,  the  specific  enthalpy  for  each  species,  /i^,  is  given  by 

h,  =  £r;dr  + hi 


(12) 


where  is  obtained  from  curve  fits  based  on  a  combination  of  experimental  data  and  theoretical 
modeling,  while  is  the  heat  of  formation  of  species  s  at  OK.  These  curve  fits  are  of  the  form 


«!,*  +  +  04, +  Os,, 


(13) 


and  the  coefficients  a,,  are  taken  from  Ref.  27.  For  the  two-temperature  model,  is  written  as 

ht  =  bir,a  +  hv^,  (^4) 


where 
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V,=  rc;^t,dr  +  h°  =  etr,»  +  ^, 

Jo  p. 


(15) 


and 


hv,a 


C‘ydr  =  ev.. 


(16) 


Here  is  the  constant  translational/rotational  portion  of  the  specific  heat  (C^ 

diatomic  species  and  2.SR^  for  monatomic  species).  The  vibrational/electronic  portion  of  the 

specific  heat  is  calculated  according  to  the  formula 


fis 

t-p.l/ 


c;{Tv)-ci^. 


(17) 


Here,  C‘(T^  means  that  the  curve  fits  in  Eq.  (13)  are  evaluated  using  Ty.  The  temperatures  T  and 
Ty  are  determined  from  and  Ey . 

The  stress  tensor  with  the  Stokes  hypothesis  is  expressed  as 


dw,  du,  ^  2  duk  - 


(18) 


and  the  components  of  the  heat  flux  vector  are  given  by 


dT 

-^Ir^ - 

dXj 


dTy 

dxj ' 


(19) 


The  individual  species  viscosities,  |X^,  are  calculated  using  the  following  curve  fits  (Ref.  27); 

fi,  =  0.1  exp(65.«  +  b4,^T  +  63,,  7"^  +  (20) 

where  T=  In  T.  The  species  conductivities  are  related  to  the  species  viscosities  by  the  relations  (Ref. 
27) 


and 


^(r,a  —  3.15fiaRf  + 


(21) 


—  P'^sj^'u.V ' 


(22) 
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where  the  binary  diffusion  coefficient  is  given  by  the  following  curve  fit: 


10.1325 


®Xp(<^4,a,T  +  +  d2.*,r7^  + 


The  mixture  viscosity  and  mixture  conductivity  are  computed  using  the  semi-empirical  methods  of 
Armaly  and  Sutton  (Refs.  28, 29). 

The  species  diffusion  velocities,  ,  are  given  by  Pick’s  law, 


The  above  equation  is  only  rigorously  valid  for  cases  in  which  the  following  conditions  hold  (Ref. 
30):  ( 1)  the  gas  is  a  binary  mixture;  (2)  the  molecular  weights  of  the  species  in  the  binary  mixture 
are  equal,  or  the  pressure  is  constant;  (3)  thermal  diffusion  is  negligible;  and  (4)  the  body  forces  per 
unit  mass  acting  on  each  species  are  equal.  If  these  conditions  are  not  satisfied,  then  the  rigorous 
computation  of  the  diffusion  velocities  involves  the  solution  of  a  matrix  equation.  However,  the 
complexity  of  this  process  has  lead  to  the  concept  of  an  effective  binary  diffusion  coefficient,  D^, 
(Ref.  30)  defined  to  be 


D, 


1-n 

V  Xn  ' 


(25) 


where  is  obtained  by  treating  ns  -  1  of  the  species  as  if  they  were  present  in  trace  amounts  and 
diffusing  through  a  background  of  a  predominate  species.  Here,  is  defined  to  be  an  average  of 
the  binary  diffusion  coefficients  T>,^  where  T>,  ,  is  the  diffusion  coefficient  obtained  for  a  mixture 
composed  of  only  species  s  and  r  (Ref.  31).  Note  that  in  order  to  strictly  enforce  mass  conservation, 
the  diffusion  velocities  of  the  predominate  species  would  have  to  be  obtained  from  the  relationship 
=  0,  j  =  1,  2,  or  3.  However,  this  expression  is  not  employed  in  the  present  code.  The 
diffusion  velocities  for  all  species  arc  computed  using  Eq.  (24).  This  expression  is  a  reasonable 
approximation  provided  that  the  molecular  weights  of  the  species  do  not  vary  widely,  which  is  the 
case  for  air. 

The  species  source  terms,  (o^,  and  vibrational/electronic  energy  source  term,  o)^  are  also 
provided  by  NEQPAK.  The  species  source  terms  are  constmcted  as  follows.  Let  the  rth  reaction 
involving  species  s  be  written  in  the  form 
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Yi  ‘'r.uAf,, 
n=l 


n$ 

^-.ni^nn 

m=l 


(26) 


in  which  M  represents  one  mole  of  species  s  and  the  and  j  are  the  stoichiometric 
coefficients.  The  rate  of  disappearance  of  a  species  s  due  to  the  rth  reaction  is 


n.t  r  nt 

£r,*  =  tK,sK  n  fl 


n=l 


TL=1 


(27) 


while  the  production  rate  of  a  species  s  due  to  the  /th  reaction  is 


UJ  ^ 

n  7/^"‘  +  Vr.sK  n  7m  ■' 


m=  1 


m=l 


(28) 


The  net  rate  of  change  of  species  s  due  to  all  reactions  is 

nr 

=  Ms  -  C,r,a), 

r=1 


(29) 


where  nr  is  the  number  of  reactions.  Like  most  aerothermochemical  models,  NEQPAK  assumes  the 
modified  Arrhenius  form  (Ref.  32)  for  chemical  reaction  rates 

H  =  ArT^^exp{-CrlT),  (30) 

where  Tis  the  translational  temperature  of  the  gas.  The  rates  defined  in  Eq.  (30)  were  obtained  from 
fits  to  experimental  data  that  were  collected  under  conditions  of  thermal  equilibrium.  However, 
when  the  reaction  rale  under  question  represents,  for  instance,  a  two-body  dissociation  reaction, 
then  the  reaction  rate  also  depends  on  the  vibrational  temperature  of  the  gas.  One  of  the  more  pop¬ 
ular  approaches  used  to  account  for  thermal  nonequilibrium  is  that  of  Park  (Ref.  33),  who  replaced 
the  translational  temperature  in  Eq.  (30)  with  a  generic  temperature  or  average  temperature,  T,.  Park 
defined  by  the  relation 

r,  =  <  Q  <  1,  (3i) 


where  a  was  chosen  to  reproduce  experimental  data. 
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The  reverse  reaction  rate  coefficient  k[  is  related  to  the  forward  rate  by 


K  =  H/k 


C 

T  » 


(32) 


where 


A  F 


(33) 


is  the  equilibrium  constant  for  the  rth  reaction  and 

n» 

5=1 

n5 

AF  = 

5=1 


(34) 


In  Eq.  (33),  p^tm  =  1  01 325  x  10*  Pa.  Also  is  the  Gibbs  free  energy  for  species  s  at  the  given 
temperature  and  a  pressure  of  \atm.  The  Gibbs  free  energies  are  computed  from  curve  fits  taken 
from  Ref.  27. 


The  vibrational  source  term,  oo^,  is  the  sura  of  two  terms.  The  first  term,  0)^  ,  models  the 
relaxation  of  vibrational  energy  under  collisions  with  heavy  particles  while  the  second  term,  , 

deals  with  the  addition  or  removal  of  energy  from  the  vibrational/electronic  energy  pool  due  to 
chemical  reactions.  The  vibration-translation  interaction  is  modeled  according  to  the  theory 
developed  by  Landau  and  Teller  (Ref.  34).  Landau  and  Teller  showed  that  if  the  vibrational  energy 
levels  of  a  molecule  are  equally  spaced  (harmonic  oscillator  approximation)  and  only  single-level 
transitions  are  allowed,  then 


(^0,5 


3=1 


(35) 


where  the  ,  indicates  the  equilibrium  value  of  the  vibrational  energy,  and  is  the  vibrational 
energy  relaxation  time.  Here,  e,,  represents  the  vibrational  energy  of  the  molecule  rather  than  the 
vibrational/electronic  energy.  Strictly  speaking,  Eq.  (35)  was  derived  for  a  harmonic  oscillator,  but 
it  is  applied  to  anharmonic  oscillators  by  using  the  polynomial  curve  fits  defined  in  Eqs.  (13)  and 
(16)  to  calculate  the  vibrational  energy  as  follows.  First,  the  electronic  excitation  energy  is 
computed  from  the  formula 

j  iipi 

Cf.s  =  M,—  53  exp(-V'(),  (36) 

^  /=i 
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where  ^i=ttlkTy,  and  the  electronic  partition  function  is 

nr/ 

1=1 


(37) 


where  nei  is  the  number  of  electronic  energy  levels  for  species  s;  is  the  degeneracy  of  the  /th 
electronic  level,  and  E,is  its  energy.  is  set  equal  to  the  difference  of  and 

The  relaxation  time  of  species  j  in  a  bath  of  species  r  is  calculated  using  the  empirical  cor¬ 
relations  of  Millikan  and  White  (Ref.  35),  where  is  defined  to  be 


pr.,,  =  exp[>l,.,7’-‘^-’  -  O.QiaAs^rfiHr  -  18.42],  (38) 

where 

/I,.,.  =O.OO116<C'0</^  (39) 

p  is  the  pressure  in  atmospheres,  and  ji,,,  is  the  reduced  mass  of  the  interacting  molecules  s  and  r. 
In  Eq.  (39),  0,  =  hvjk,  where  v,  is  the  characteristic  vibrational  frequency  of  species  s\  h  is  the 
Planck  constant,  and  k  is  the  Boltzmann  constant.  The  relaxation  time  employed  in  Eq.  (35)  is 
computed  from  the  'Cg.^  according  to  the  formula 


where 


US 


r=l 


(40) 


(41) 


is  Park’s  high-temperature  correction  (Ref.  36)  and 

C  =  =  2.89S677  x  10“^ 


in  SI  units.  The  value  =  10  “  ra*  for  the  effective  cross  section  for  vibrational  relaxation  is  that 
suggested  by  Gnoffo  et  al.  (Ref.  26).  The  second  term  in  taj,  models  the  effect  that  chemical 
reactions  have  on  the  vibrational  energy.  Clearly,  the  dissociation  of  a  molecule  removes 
vibrational/electronic  energy,  while  recombination  of  a  molecule  must  add  vibrational/electronic 
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energy  to  the  species  pool.  In  this  work  it  is  assumed  that  nonpreferential  dissociation  occurs.  Under 
this  assumption,  takes  the  form 


ns 

(Jjy  = 

1 


(42) 


3.0  NUMERICAL  FORMULATION 


3.1  THREE-DIMENSIONAL  FINITE-VOLUME  FORMULATION 


The  following  derivation  is  for  the  two-temperatuie  model.  However,  the  equations  and 
numerical  method  are  easily  reduced  to  the  one-temperature  model.  The  integral  formulation  of  a 
system  of  conservation  equations  over  a  finite  region  of  space  takes  the  general  form  (Ref.  37) 


£ 

dt 


+  IS/-^ 


(43) 


where  for  flows  in  thermo-chemical  nonequilibrium 


/  Pi  'i 

Pua 

Ev 

;  n  = 

OJy 

pu 

0 

pv 

0 

pw 

0 

E  J 

1  0 

(44) 


and  the  numerical  flux  vector  is  written  as  the  sum  of  inviscid  and  viscous  contributions; 

=  (Fl  -h  /’v)tr  +  (fVl  +  Cwy)iy  +  (^Tl  -|-  Hv)ig,  (45) 

Q  and  £2  are  the  vectors  of  conserved  quantities  and  source  terms,  respectively.  F,  G,  and  H  are  the 
inviscid  flux  vectors.  They  are  equal  to  the  flux  of  Q  per  unit  area  per  unit  time  in  the  Cartesian 
coordinates  x,  y,  and  z,  respectively,  with  unit  vectors  lx,  ly,  and4 .  For  the  present,  we  will  ignore 
the  viscous  flux  vectors  F„  G„  and  Hy.  The  inviscid  flux  vectors  may  be  written  as 
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Pi « 

^  plV  ^ 

i  p\vi 

Pn^U 

PnsV 

PnsW 

Ev  u 

;  fv,  = 

Ev  tJ 

;  fii  = 

Ev  w 

pu^  +  p 

ptiv 

puw 

pvu 

pv^  +  p 

pvw 

/Mjm 

pwv 

pw^  +  p 

iE  +  p)u  ) 

{(E  +  p)v  ) 

^  {E  +  p)w  J 

(46) 


In  order  to  solve  Eq.  (43)  numerically,  the  spatial  domain  of  interest  is  broken  into  a  set  of 
disjoint  control  cells  with  volumes,  V,  and  cell  face  areas.  Oj.  The  flux  into  or  out  of  a  finite^ell 
volume,  V,  across  a  cell  face  Oj  is  obtained  by  a,- .  The  equations  arc  written  in  terms  of  the 
fluxes  in  the  computational  directions  F,  G,  and  H.  These  directions  correspond  to  the  integer  indi¬ 
ces  j,  k,  and  /. 


An  implicit  fuiite-volume  discretization  of  Eq.  (43)  is 


At 


(47) 


Throughout  this  report,  the  notation  of  Gnoffo  (Ref.  21)  has  been  adopted  where  capital  letters 
J,  K,  L  indicate  cell  centers  and  j,  k,  I  indicate  cell  faces.  Also,  unlike  most  finite— volume  schemes, 
NED  ANA  stores  the  conserved  quantities  at  the  cell  vertices.  The  difference  between  this  technique 
and  traditional  finite-volume  techniques  is  largely  transparent  (see  Appendix  B).  Storing  the 
conserved  variables  directly  at  mesh  points  allows  for  interaction  with  finite-difference  algorithms 
without  altering  the  computational  grids  or  data. 


The  numerical  flux  function  is  based  on  the  flux-limited,  artificial  dissipation  model  of  Jame¬ 
son  (Ref.  38),  and  Yoon  and  Kwak  (Ref.  39).  This  model  combines  the  simplicity  of  artiricial 
dissipation  schemes  with  the  total  variation  diminishing,  (TVD),  property.  A  discussion  of  the 
relevant  numerical  issues  involved  in  the  selection  of  this  particular  flux  model,  as  well  as  a  com¬ 
plete  description  of  the  numerical  flux  function,  is  included  in  Appendix  C.  The  flux  at  the  cell  face 
y  +  1  is  calculated  as  the  average  of  the  fluxes  from  cell  J  and  7  -t- 1 .  Ihis  averaging  requires  the  flux 
function  to  be  evaluated  twice  and  then  averaged.  An  alternative  approach  would  be  to  average  the 
conserved  variables  to  obtain  values  at  the  cell  face  and  then  perform  the  flux  evaluation.  The 
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former  method  was  chosen  for  its  flexibility  in  writing  finite-difference-type  schemes.  The  numer¬ 
ical  flux  in  the  ^  direction  at  the  j  +  I  face  has  the  form 


r’l+i 


-  2.0A^Q'}X\l 


(48) 


where  the  cell  face  directed  area  is  defined  as 

Cf^  =  +  (T|i„ 


(49) 


and 


«+l 


Here,  4*  is  defined  for  vectors  P  and  Q  to  have  components: 

4^,1/*, (?)  =  5[sigii(/",)  +  sigu((5,)]iiiin(|P,-|,|Q,|). 


The  scalar  dissipation  coefficient  is  defined  by 

C^J+1.A',L  =  0.2.5[A^-j+i,a'.L  +  j.k',£,][k2  +  K4 J.K'.l], 

and  the  spectral  radius  of  the  flux  Jacobian  is  defined  at  the  cell  center  by 


J,K,L 


(T^ 


J.h\L 


(50) 


(51) 


(52) 


(53) 


and  the  cell  area  in  the  ^  direction  evaluated  at  the  cell  center  is 

— 

f'iJ.K.L  -  +  ^ihK,L)- 


(54) 


The  frozen  speed  of  sound  is  defined  as 


(55) 
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where  the  frozen  ratio  of  specific  heats  is  expressed  as 


=  (1  +  13), 


(56) 


where  P  =  3p/3E  (see  Appendix  D). 

The  dissipation  Is  a  function  of  the  pressure  sensor 

(57) 


where 


P./+l,A'.t<  -  'i-^PJ,K,L  +  pJ-\,K.L 
PJ+\M.L  +  '^0PJ,K,L  +  PJ-1,A,L 


(58) 


The  numerica)  fluxes  C,  and  H  are  calculated  similarly. 

3.2  QUASI-ONE-DIMENSIONAL  FORMULATION 

The  one-dimensional  form  of  the  conservation  equations  with  area  change  is  derived  in 
Appendix  E.  This  form  of  the  equations  is  solved  when  performing  numerical  studies  for  shock  tube 
and  nozzle  problems.  The  finite-difference  form  of  the  equations  is  presented  below: 


/ 

■  Pi 

( 

\ 

■  Ay-'wi  ■ 

AJ-^ 

Pna 

+ 

AJ-^ 

Ev 

EyV^x 

AJ~^u}v 

pu 

AxJ~^p 

E 

} 

t 

\ 

i 

0 

Note  that  the  equation  has  been  transformed  to  curvilinear  coordinates.  The  Jacobian  of  the 
transformation  is  7  =  d^/dx.  For  constant  area  problems  such  as  shock  tubes,  A  is  set  equal  to  I . 

3.3  VISCOUS  TERMS 

In  this  section  the  viscous  contributions  to  the  flux  function  will  be  discussed.  These 
contributions  include  shear  stress,  heat  conduction,  and  diffusion  terms.  The  viscous  Cartesian  flux 
function  for  the  x  direction  from  Eq.  (4S  )  has  the  form 
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^  Pi  «‘{ 

Pn, 

HaLlPtliVs'ui  +  9Vi 


+  qr  +  qvx 

\  +T:ctU  +  rj:yV  +  Tj.gVr  } 


(60) 


Each  component  of  Eq.  (60)  may  be  fonnulated  in  the  context  of  a  finite-volume  scheme.  A 
thin-layer  approximation  also  is  employed  in  the  direction  normal  to  the  body  surface,  Recalling 
that  f  -Tj  ,  the  viscous  flux  contribution  for  species  continuity  becomes, 


where 


The  viscous  flux  contribution  to  the  vibrational — electronic  energy  equation  can  be  expressed  as 


-I>s 

dc 


J  J.K.t 


(63) 


where 


/  aTv\ 

■"'t  =  “('“"ac'Ar.w'  («4) 

The  viscous  flux  contribution  to  the  x-momentum  equation  in  finite-volume  form  can  be  written  as 


=  -p 


'  na+U 


du  ^  \  du  ^  . 


J,K,t 


(65) 
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In  a  similar  manner,  the  viscous  flux  contribution  for  all  of  the  momentum  equations  can  be  written 
compactly  as 


ti»+4 


+ 

+ 

+ 


^2  G  \ 

P-B-zCy 


(66) 


where 


and 


(67) 


(68) 


The  viscous  flux  contribution  to  the  total  energy  equation  now  can  be  easily  written  using  the 
previously  developed  expressions.  The  equation  becomes 


^  L  dXa 

tftrC  +  9V  (  +  2^  Paf^a 


\  (  du  dv  \  1 


where 


+  /i[/72  (uG  +  uG  + 

3A  CHIMERA  DOMAIN  DECOMPOSITION 


(69) 


(70) 


NEDANA  has  incorporated  the  chimera  domain  decomposition  procedure  developed  by 
Benek,  Steger,  and  their  co-woikers  (Refs.  7, 40-43).  The  chimera  scheme  was  developed  to  allow 
a  system  of  simple  grids  with  simple  topologies  to  model  complex  aerodynamic  configurations  or 
bodies  in  relative  motion. 


The  general  concept  behind  chimera  is  illustrated  in  Fig.  1,  which  depicts  two  independently 
generated  meshes  representing  a  flapped  airfoil.  The  flap  mesh  is  embedded  within  the  airfoil  mesh. 
Clearly,  the  flap  mesh  outer  boundary  can  receive  flow-field  information  interpolated  from 
appropriate  mesh  points  of  the  airfoil  mesh.  However,  a  reverse  process  must  occur  as  well;  the 
airfoil  mesh  must  receive  flow-field  information  from  the  flap  mesh.  This  transfer  is  achieved  as 
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follows.  A  hole  must  be  created  within  the  airfoil  mesh  to  remove  the  points  within  the  airfoil  mesh 
that  are  interior  to  the  flap.  Hie  hole  boundary  points  of  the  airfoil  mesh  can  be  updated  by 
interpolation  from  the  flap  mesh.  In  general,  any  mesh  can  receive  information  from  other  meshes 
through  both  outer  boundary  and  hole  boundary  points. 


Flap  Mesh 

Figure  1.  Mesh-to>mesh  communication. 

The  interpolation  process  is  further  illustrated  in  Fig.  2,  which  depicts  a  portion  of  the  overlap 
region  between  the  airfoil  and  flap  meshes.  Airfoil  mesh  points  inside  the  hole  region  surrounding 
the  flap  are  blanked  out  of  the  computational  domain  of  the  airfoil  mesh.  In  chimera  terminology 
they  are  hole  points.  The  hole  region  is  defined  by  a  hole  creation  boundary  within  the  flap  mesh. 
The  points  in  the  airfoil  mesh  surrounding  the  blanked  points  are  hole  boundary  points,  and  they 
receive  flow-field  information  interpolated  from  mesh  points  within  the  fl^  mesh. 
Correspondingly,  points  on  the  the  outer  boundary  of  the  flap  mesh  receive  flow-field  information 
interpolated  fiom  mesh  points  within  the  airfoil  mesh. 

Application  of  the  chimera  scheme  requires  two  steps:  (1)  a  description  of  how  each  mesh  is 
to  communicate  flow-field  information  to  other  meshes,  and  (2)  execution  of  the  flow  solver  (in 
this  case,  NEDANA)  using  the  information  generated  in  step  1.  Step  1  is  performed  by  PEGSUS 
(Ref.  44)  developed  at  the  AEDC.  The  processes  accomplished  by  PEGSUS  include  establishing 
which  boundary  points  in  a  mesh  will  be  updated  by  interpolated  flow  variables  from  other  meshes. 
Also,  PEGSUS  calculates  the  required  interpolation  coefficients  for  donor  mesh  points  that  provide 
interpolated  information  for  recipient  boundary  points  in  another  mesh.  More  information  on  the 
chimera  scheme  and  the  PEGSUS  code  can  be  found  in  the  references. 
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Airfoil  Mesh 


interpolated  Boundary  Point 
In  Flap  Mesh 


Interpolated  Boundary 
Point  in  Airfoil  Mesh 


Hole  Creation 
Boundary  in 
Flap  Mesh 
Hole  Boundary  In 
Airfoil  Mesh 


Hole  in  Airfoil 
Mesh 


Flap  Mesh 


Figure  2.  Overlap  region  between  meshes. 


3.5  SOLUTION  PROCEDURE 


Equation  (47)  is  a  set  of  nonlinear  algebraic  equations  which  is  solved  at  each  time  step.  A 
variety  of  techniques  exists  to  solve  such  equations.  An  encyclopedic  review  of  such  techniques 
may  be  found  in  Ref.  4S.  Throughout  this  section,  the  nomenclature  of  Ref.  4S  is  used. 


Newton's  method  is  the  classical  procedure  for  solving  nonlinear  algebraic  equations.  For 
this,  first  define  the  residual  vector,  /?"(g),  at  time  level  n,  as  a  function  of  a  vector,  g,  so  that  at  a 
point  J,  K,  Ly  R”  has  components 


iQj,K,L-Qj,K.L) 

=  ^ - 

+  {f^'j.k+i.L(Q}  -  f’JxUQ)) 

+  (Hj.K.t+iiQ)  -  ffj,Kj{Q)] 


The  solution  to  Eq.  (47)  at  time  level «  +  1  is  the  vector,  g"*^that  satisfies  the  residual  equation. 


R’‘(g’‘+')  =  o. 


(71) 
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An  approximate  Newton’s  method  is  used  to  solve  this  equation  iteratively.  The  initial  iteration  is 
always  ^  sufficient  number  of  iterations,  0**^  is  set  equal  to  the  last  iteration.  To 

sinnplify  notation  in  the  following,  the  dependence  on  n  is  suppressed.  Let  denote  the 
iterative  approximation  to  the  solution  of  Eq.  (71).  Also,  let  denote  the  residual  vector 
evaluated  at  Similarly,  let  the  evaluation  of  the  fluxes  and  the  sources  at  be  denoted  by 
the  superscript,  (m). 


To  motivate  the  development  of  the  iterative  techniques  used  in  the  NEDANA  code,  consider  the 
use  of  a  pure  Newton's  method  approach  to  solving  Eq.  (71).  First,  the  following  linear  system  of 
equations  is  solved; 


(72) 


or 


g(m+l) 


(p(m) 


R<"‘^ 


(73) 


Next,  the  norm,  -  Q”>\\,  is  checked  to  see  if  it  is  less  than  some  prescribed  tolerance.  If  it  is, 

then  is  accepted  as  the  solution  to  Eq.  (47)  at  time  level «  + 1 .  If  I  is  greater  than 

the  prescribed  tolerance,  then  m  is  incremented  and  Eq.  (73)  is  solved  again.  This  process  is  repeat¬ 
ed  until  the  norm  is  less  than  the  prescribed  tolerance  or  until  a  fixed  number  of  iterations  has  been 
performed.  When  solving  problems  in  two  and  three  dimensions,  Newton’s  method  is 
computationally  intensive  since  the  Jacobian,  is  a  very  large,  sparse  matrix.  This  makes 

direct  inversion  of  the  Jacobian  impractical.  Therefore,  as  described  in  detail  below,  the  NEDANA 
code  uses  an  alternative  to  Newton’s  method  to  solve  Eq.  (71).  This  method  is  referred  to  as  a  mod¬ 
ified  one-step,  odd/even,  Jacobi-Gauss-Seidel-Newton  (JGSN)  scheme.  In  the  past,  this  iterative 
method  was  known  as  a  variant  of  the  LIM  scheme. 

To  motivate  the  use  of  the  JGSN  method,  first  examine  the  Gauss-Seidel-Newton  (GSN)  and 
the  Jacobi-Newton  (JN)  schemes.  If  a  one-step  GSN  iteration  were  employed  to  solve  Eq.  (71),  a 
Gauss-Seidel  method  would  be  applied  as  a  primary  iteration  and  Newton’s  method  would  be 
applied  as  a  secondary  iteration.  Specifically,  in  the  primary  Gauss-Seidel  iteration,  for  every  grid 
point,  J,  K,  L,  the  local  residual  equation. 


(74) 
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J  I  \  ^ 

is  solved  for  Qj  ,i  i_  .  Here,  m  indicates  that  the  components  of  g‘'»*‘>are  used  in  Eq.  (74)  if  they 
are  available,  and  otherwise,  the  components  of  are  used.  Then,  in  the  secondary  iteration, 
Newton’s  method  is  used  to  solve  Eq.  (74)  according  to 

^J.K,L  Qj,K.L 

In  a  one-step  version  of  this  scheme,  only  one  iteration  of  Eq.  (75)  is  implemented  before  proceed¬ 
ing  to  the  next  grid  location.  When  all  components  of  are  calculated,  a  sweep  is  said  to  be 
completed. 

If  a  one-step  JN  iteration  were  employed  to  solve  Eq.  (71),  a  Jacobi  method  would  be  applied 
as  a  primary  iteration  and  Newton’s  method  would  be  applied  as  a  secondary  iteration.  Specifically, 
in  the  primary  Jacobi  iteration,  for  every  grid  point,  J,  K,  L,  the  local  residual  equation. 


p(TFl*) 

(75) 


E'J, A", Q 


{,«) 

J,A',L-1 


'1  -  n 


(76) 


is  solved  for  •  Note  that  here  the  fixed  arguments  of  the  local  residual  are  written  only  in 

terms  of  the  components  of  Then,  in  the  secondary  iteration,  Newton’s  method  is  used  to  solve 
Eq.  (76)  according  to 


9<?Kx 


1  -1 


j»{m) 


(77) 


Both  methods  described  above  involve  only  the  inversion  of  an  x  matrix,  where  nq  is 
the  number  of  conserved  variables.  However,  the  GSN  method  is  recursive  and  is  not  well  suited 
for  implementation  on  machines  that  employ  vector  processors  (Ref.  46).  In  addition,  Tramel  (Ref. 
16)  has  shown  that  a  JN  iteration  is  unstable  for  CFL  numbers  greater  than  one.  Therefore,  an  odd/ 
even  JN  iteration  is  proposed  in  Ref.  16  and  is  the  basis  for  the  solution  technique  used  in 
NEDANA. 


In  NEDANA,  an  odd/even  JN  iteration  is  applied  along  lines  of  constant  J,  K,  while  a  GSN 
iteration  is  ^plied  in  the  other  computational  directions.  As  indicated  above,  this  procedure  is 
referred  to  as  a  Jacobi-Gauss-Seidel-Newton  (JGSN)  method.  The  odd/even  JGSN  iteration  is 
implemented  in  NEDANA  as  follows.  First,  the  following  equation  is  solved  for  every  point  along 
a  line  of  constant  J,  K,  with  L  =  L^,  an  even  value. 
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Here, 


/)(«»+*)  _ 

Vj,K,Lr  -  '‘iJ.K.U  ~ 


(m) 


n  -1 


D  (’«*).'■ 


(78) 


l.K,L, 


~  “i.K.LrJ 


+ 

+ 


r{«n  >1 


+  Vj,A-.L, 


At 


-  V 


(m*) 


(79) 


Again,  m*  indicates  that  the  components  of  are  used  if  they  are  available,  and  otherwise,  the 
components  of  are  used.  In  addition,  the  superscript,  e,  in  emphasizes  that  H  is 

evaluated  entirely  in  terms  of  Qf"*K  as  opposed  to  being  evaluated  as  described  below  for  the  odd 
Jacobi  iteration.  Note  that  the  Jacobian  matrix,  o''®’’  course  of  a 

time  step.  Next,  the  following  equation  is  solved  for  every  point  along  a  line  of  constant  J,  K,  with 
L  =  L^,  an  odd  value,  ^ 

dR  I 


Vj.A'.Lo 


^■J,K,Lo  ~ 


^JJ<,Lo 


-t 


\  ‘^J+^,K.U 


R 


J,K,U- 


-  F\ 


r(™*) 


J.A'Xo 


(80) 


Here, 


+ 

+ 


f  (»'‘),o 


“  'V.fc,Z,c 


) 


~  Q'j,k,u  ) 


+  yj,KX 


At 


_  V 


i1 


(m) 


(81) 


Recall  that  ( 1  -  depends  on  {2y,#:,t.2*0y.r.i,^i«l2j.iu-»G/,x.i*i»G/.x,Lta}-Iuthiscase,  the 

superscript,  o,  in  emphasizes  diat  is  evaluated  in  terms  of  when  L  is  odd,  and 

in  terms  of  ^ 
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Once  all  points  on  the  line  of  constant  J  and  K  have  been  updated,  the  iteration  proceeds  to 
the  next  line.  When  all  components  of  are  calculated,  a  sweep  is  said  to  be  completed. 

Typically,  the  total  number  of  sweeps,  ,  is  set  to  two  or  four,  and  0^“'^  is  accepted  as  the  solution 
to  Eq.  (47)  at  time  level  n  +  1 . 

Computational  experience  and  theoretical  analysis  (Refs,  12-16)  have  shown  that  the  use  of 
the  exact  Jacobian,  dR  JdQ ^ ,  leads  to  an  unstable  iteration  scheme  when  large  CFL  numbers 
are  used.  In  the  past,  several  ad  hoc  formulations  have  been  advocated  on  how  to  approximate  the 
Jacobian  with  another  matrix  However,  viewing  the  LIM  scheme  as  a  Diagonalized 
Approximate  Newton’s  Algorithm  (DANA),  convergence  of  the  method  is  guaranteed  if  the 
Jacobian  of  the  iteration  process  has  spectral  radius  less  than  or  equal  to  one  (Ref.  45).  The 
importance  of  this  new  perspective  of  the  LIM  scheme  is  that  a  well-founded  criterion  is  available 
to  estimate  under  what  conditions  the  scheme  is  convergent.  This  criterion  allows  the  previous  ad 
hoc  formulations  of  to  be  readily  checked. 

In  previous  applications  of  the  LIM  scheme,  the  matrix  B  was  taken  to  be 


where 


-I- |o.r25  -f-  2.0A;ej  -I- 

-|-1.5(cfj+,  + 

-I- [o.125{A„k+i  +  2.0A„/i-  +  A„k-_i) 
+  1  .i5  (c, I  + 

+  [0.12,'3(A(;£,+  ]  -f-  2.0  A,;/,  +  A{;i,_i) 

-|-l..'5(Cf(+i  -I-  £(;{}]  +  bv. 


(83) 


and 


(84) 
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where  is  a  relaxation  parameter  optimized  for  stability  and  convergence.  However,  in  the  present 

work  the  matrix  B  becomes 


where 


Bi 

0 


0 

B2 


J,K,L 


Bi  j./f.L 


(ii»+we)x(w*+tir) 


da 

dQ 


J,K,L 


(85) 


(86) 


and 


^2J,K,L  =  l4x4. 

Here  ns  and  ne  are  the  number  of  nonequilibrium  species  and  energies,  respectively.  Note,  this 
partitioning  of  B  ignores  the  derivatives  of  w  with  respect  to  pu,  pv,  pw,  and  E.  However,  it  has  been 
observed  that  this  simplirication  does  not  affect  the  stability  or  accuracy  of  the  scheme.  In  Appendix 
F,  the  complete  derivation  of  the  nonequilibrium  source  Jacobian,  dCi/  dQ,  is  included. 

Thus,  the  nonequilibrium  variables  (pp...,  p^^,  E^^aie  solved  for  by  a  matrix  inversion,  and 
the  variables  (pw,  pv,  pw,  are  solved  for  by  a  scalar  inversion.  Speciflcally,  at  the  beginning  of 
a  given  time  step,  an  LU  decomposition  of  B/j/^ j^is  performed  for  every  point.  The  decomposition 
is  saved  and  not  updated  during  a  sweep.  In  fact,  in  the  absence  of  significant  temporal  variations, 
the  decomposition  need  not  be  updated  at  every  time  step.  Once  the  higher  order  work  of  an  LU 
decomposition  is  performed,  sweepwise  updates  of  the  nonequilibrium  variables  can  be  computed 
with  the  lower  order  work  associated  with  a  simple  back  substitution. 

4.0  RESULTS  OF  ONE-DIMENSIONAL  COMPUTATIONS 

Extensive  one-dimensional  numerical  tests  were  performed  with  NEDANA  to  compare  the 
results  of  the  present  numerical  technique  with  the  results  of  existing  flow  solvers.  Three  cases  were 
chosen: 

1.  Shock  tube  problem. 

2.  Supersonic  duct  flow  with  area  change. 

3.  Supersonic  duct  flow  with  area  change  and  normal  shock. 
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For  at]  results  presented  in  this  section,  the  one-temperature  model  was  assumed.  One-dimensional 
computations  using  the  two-temperature  capabilities  of  NEDANA  may  be  found  in  Ref.  18.  These 
three  cases  are  used  to  test  the  speed,  accuracy,  and  stability  of  the  present  algorithm  as  compared 
with  state-of-the-ait  flow  solvers  from  other  researchers. 

In  the  comparisons  presented  in  Figs.  3-8,  variables  are  normalized  by  reference  quantities.  In 
all  cases,  the  reference  quantities  are  denoted  by  an  superscript.  This  should  not  be  confused  with 
the  common  usage  of  the  *  superscript  to  indicate  sonic  conditions.  For  the  shock  tube  problem,  these 
variables  are  normalized  with  the  equilibrium  properties  of  the  high-pressure  side  of  the  di^hragro. 
The  variables  are  normalized  by  the  properties  at  the  supersonic  inlet  for  the  two  duct  problems. 

4.1  SHOCK  TUBE 

Shock  tubes  are  long  constant-area  devices  that  are  initially  pressurized  with  two  gases 
separated  by  a  diaphragm.  The  gases  are  maintained  in  states  of  equilibrium  that  are  very  much 
different  from  one  another.  For  instance,  in  the  tube  shown  in  Fig.  3,  the  left  side  of  the  tube  is 
maintained  at  a  pressure  of  100  atm  and  a  temperature  of  9000  K.  The  right  side  is  at  300  K  and  the 
pressure  is  1  atm.  The  flow  in  the  shock  tube  can  be  modeled  as  a  one-dimensional  inviscid  wave 
propagation  problem.  At  r  =  0,  ^e  diaphragm  is  burst,  and  a  strong  compression  wave  (shock) 
moves  to  the  right  with  velocity  C.  As  the  shock  wave  propagates  to  the  right,  an  expansion  wave 
will  propagate  to  the  left  as  the  mixture  responds  to  the  the  adjacent  fluid  moving  to  the  right.  Also, 
a  contact  surface  separating  the  two  gases  will  follow  the  shock  wave  as  the  shock  wave  moves  to 
the  right  compressing  the  fluid  it  encounters.  The  contact  surface  travels  with  velocity  V  where  V 
<  C.  Figure  4  details  the  various  structures  present  as  they  advance  in  time.  Notice  that  the  contact 
surface  is  maintained,  since  there  is  no  mechanism  to  mix  the  gases.  In  the  following  shock-tube 
comparisons,  relative  to  the  distance  the  shock  wave  has  traveled,  time  has  not  advanced  far  at  all. 
Consequently,  the  contact  surface  is  close  behind  (slightly  to  the  left  of)  the  shock  wave.  The  double 
step  is  most  noticeable  in  plots  of  the  temperature.  This  is  because  the  shock  wave  increases  the 
temperature  in  advance  of  the  contact  surface.  Most  of  the  temperature  increase  then  occurs  across 
the  contact  surface  for  this  particular  case.  Note  that  the  contact  surface  does  not  compress  and  heat 
the  gas  like  the  shock  wave.  The  changes  across  the  contact  surface  are  embedded  in  the  fluid  due 
to  the  initial  differences  in  the  states  of  the  gas  that  existed  across  the  diaphragm.  Also,  keep  in  mind 
that  there  is  no  velocity  or  pressure  jump  across  the  contact  surface.  Therefore,  by  observing  the 
velocity  and  pressure  curves,  one  can  determine  how  well  an  algorithm  is  capmring  shocks.  Also, 
by  observing  the  temperature  and  density  plots,  one  can  see  both  the  shock  and  contact  surface  most 
clearly. 

The  shock-tube  problem  is  ideal  as  a  first  case,  because  it  is  inherently  unsteady.  Spurious  pro¬ 
duction  of  various  constituents  of  the  gas  would  indicate  that  the  solver  was  unable  to  maintain  time 
accuracy.  To  make  a  comparison,  the  conditions  of  the  shock  tube  case  performed  by  Shuen,  Liou, 
and  van  Leer  (Ref.  47)  were  used  as  starting  conditions. 
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For  X  e  [0.0,  5.0]  rru  =  (100  aim,  9000  K,  0.0  mfsec). 

Forx  €  [5.0,  10.0]  cm  [ph.Tr,ur)  =  (1  «lm,  300  A",  0.0  mjsec). 


State  L 


Initial  State  at  f  =  0 

Initial  Position 


Expansion  Contact  Shock 

Fan  Discontinuity 

Flow  State  at />0 

Figures.  Shock  tube  sdiematic. 
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t 
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Temperature  Mixture  Pressure  Velocity 
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c.  Temperature  (normalized  by  T*) 
Figures.  Shock  tube  comparisons. 


33 


Mole  Fraction  Mole  Fraction  Mixture  Oensi 


AEDC-TR-94-18 


d.  Mixture  density  (normalized  by  p*) 


e.  Atomic  nitrogen 


f.  Atomic  oxygen 
Figures.  Continued. 
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The  equilibrium  air  densities  were  computed  using  the  curve  fits  of  Prabhu  and  Erikson  (Ref. 
48)  to  be  pL  =  2.627  kgfn^  for  the  high-pressure  side  and  pR  =  1.1733  for  the  low-pressure 
side.  In  the  following  comparisons,  all  computations  were  performed  on  a  grid  consisting  of  2(X) 
equally  spaced  nodes  distributed  over  the  interval  [0.0, 10.0]  cm.  The  comparison  presented  in  Fig. 
5  is  made  with  Shuen’s  MUSCL-interpolated  Roe  upwind  scheme  (Ref.  47)  and  with  Molvik’s  flow 
solver  (Ref.  49).  Both  Molvik's  flow  solver  and  the  present  flow  solver  were  run  for  2S0  time  steps 
at  a  CFL  number  of  O.S,  which  gave  an  absolute  time  of  0. 1 67  x  10'^  sec.  This  nun^r  of  steps  was 
chosen  so  as  to  match  the  shock  position  of  Shuen  et  al.,  who  supplied  no  temporal  data.  Boundary 
conditions  for  this  problem  are  unimportant,  as  the  flow  at  both  boundaries  remains  undisturbed. 
The  results  from  NEDANA  and  the  results  from  the  Shuen  et  al.  flow  solver  were  computed  using 
Dunn  and  Kang’s  (Ref.  50)  reaction  model  (ionization  reactions  were  not  included).  The  results 
from  Molvik’s  flow  solver  were  obtained  using  Blottner’s  (Ref.  51)  reaction  model  and  were 
obtained  at  the  AEDC  by  the  authors.  The  results  attributed  to  Shuen  et  al.  were  manually  digitized 
from  plots  taken  from  the  cited  article.  This  accounts  for  waviness  present  in  their  results. 

Similar  results  for  the  gas  constituents  and  flow  properties  are  presented  in  Fig.  5  for  all  three 
flow  solvers.  The  pressure  comparison,  Fig.  5b,  demonstrates  that  NEDANA  produces  shock  fronts 
that  are  as  sharp  as  those  predicted  by  the  flow  solvers  of  both  Molvik  and  Shuen  et  al.  Figure  5c, 
however,  shows  that  NEDANA  does  not  reproduce  the  contact  surface  as  well  as  the  other  codes.  ‘ 
This  is  because  both  the  Molvik  and  Shuen  et  al.  flow  solvers  are  based  on  approximate  Riemann 
solvers,  which  force  the  net  flux  to  be  modeled  exactly  across  discontinuities  in  one  dimension.  In 
multi-dimensional  flows  in  which  the  grid  is  not  aligned  with  the  discontinuity,  much  less  favorable 
results  are  obtained  using  Roe’s  Riemann  solver.  It  is  encouraging  that  the  present  solver  so  closely 
follows  the  Molvik  solution,  as  it  is  a  well  validated  and  widely  used  flow  solver.  The  agreement 
in  the  expansion  region  between  the  NEDANA  and  the  Molvik  results  on  the  one  hand,  and  both 
disagreeing  with  the  results  of  Shuen  etal.  on  the  other  hand,  would  indicate  that  the  solver  of  Shuen 
et  al.  does  not  handle  the  expansion  correctly.  Overexpansions  are  typical  with  flow  solvers  using 
higher  order  upwind  schemes;  but,  since  Molvik  also  uses  a  high  order  upwind  method,  it  is  not 
clear  what  mechanism  is  driving  this  disagreement. 

Mole  fractions  for  the  species  N,  O,  N.^,  and  NO  as  computed  by  the  three  flow  solvers 
are  plotted  together  in  Figs.  5e-i.  Again,  the  three  flow  solvers  produce  similar  results,  with  the 
codes  of  Molvik  and  Shuen  et  al.  producing  sharper  discontinuities  and  the  results  predicted  by 
the  Molvik  and  NEDANA  flow  solvers  agreeing  more  closely,  in  general.  Notice  that  the  relative 
amount  of  the  species  changes  greatly  across  the  contact  surface,  but  little  across  the  shock.  In 
particular,  the  mixture  on  the  high-pressure  side  has  little  molecular  oxygen.  Fig,  5h.  However, 


1.  Additional  calculations  have  been  made  with  NEDANA  using  a  doubly  finer  mesh.  The  results  show  a 
much  sharper  contact  surface,  and  the  density  comparison  is  nearly  exact. 
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because  of  the  temperature  drop  across  the  contact  surface,  the  amount  of  molecular  oxygen  rises 
and  remains  nearly  constant  through  the  shock.  This  behavior  illustrates  the  finite-rate  nature  of  the 
chemical  process.  The  motion  of  the  shock  is  very  rapid;  thus,  the  chemical  changes  lag  behind  the 
shock.  Consequently,  changes  in  the  species  concentrations  resulting  from  the  shock  are  not 
observed  until  the  changes  resulting  from  the  contact  discontinuity  begin.  Changes  across  the 
contact  surface  then  become  dominant  and  obscure  the  much  smaller  changes  initiated  by  the  shock 
wave.  By  definition,  the  molar  fractions  of  the  constituents  of  the  mixture  must  change  discontinu- 
ously  across  the  contact  surface.  This  is  not  captured  in  the  numerical  simulation.  Numerical  smear¬ 
ing  of  the  contact  surface  causes  nonphysical  intermediate  temperatures  that  are  conducive  to  the 
production  of  species  not  actually  present  at  the  contact  surface.  In  particular,  notice  that  all  the  flow 
solvers  smear  the  contact  surface  and  produce  spurious  overshoots  in  NO  and  O,  This  overshoot  is 
related  to  the  fact  that  the  equilibrium  mole  fractions  of  NO  and  O  as  a  function  of  temperature  for 
a  given  pressure  are  non-monotone.  This  spurious  production  itself  causes  further  smearing  because 
there  also  is  a  lag  in  this  nonphysical  production  of  species.  The  numerical  and  nonequilibrium 
issues  thus  become  entangled.  By  turning  off  the  chemistry,  NEDANA  produces  solutions  with  less 
smearing  of  the  species  mole  fraction  profiles  across  the  contact  surface;  but  with  the  chemistry 
enabled,  the  smearing  is  amplified. 

4.2  SUPERSONIC  DUCT 

Supersonic  ducts  with  area  change.  Fig.  6,  are  good  model  problems  for  testing  a  time¬ 
marching  code’s  ability  to  achieve  an  accurate  steady-state  solution.  Two  cases  were  chosen  for 
comparison;  ( 1)  a  supersonic  duct  flow  that  remains  supersonic  to  the  exit;  and  (2)  a  supersonic  duct 
flow  that  has  had  its  exit  pressure  raised  until  a  shock  is  standing  approximately  half-way  down  the 
duct,  with  subsonic  flow  in  the  remaining  portion  of  the  duct  from  the  shock  to  the  exit. 

Results  of  the  two  nozzle  flow  problems  are  presented  and  compared  with  the  results 
produced  by  Molvik's  flow  solver  for  the  case  with  no  shock.  Both  the  flow  solvers  use  the  same 
reaction  models  used  for  the  shock  tube  example.  However,  for  the  supersonic  duct  test  cases, 
ionization  reactions  are  included  in  the  models.  Results  are  also  presented  for  the  two  test  cases  of 
equilibrium  air  solutions  produced  using  curve  fits  taken  from  Ref.  48.  The  area  distribution  of  the 
duct  is  taken  as 


A(x)  =  .5.5  +  4.5taiili(0.7x  -  8.5);  r.  G  [0.0,  10.0]  m, 

which  produces  a  nozzle  with  an  area  ratio  of  10.0.  Both  grids  consisted  of  101  uniformly 
distributed  points.  The  supersonic  inflow  conditions  for  both  cases  were: 


7  "  =  6000  K 
=  10  atm, 

while  the  frozen  Mach  number  at  the  inflow  was  2.5. 
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In  both  of  the  supersonic  duct  test  cases,  convergence  to  steady  state  was  assnined  when  the 
unsteady  residual  was  driven  to  machine  zero  everywhere. 


Figure  6.  Duct  schematic. 

4.3  SUPERSONIC  DUCT  WITH  AREA  CHANGE 

As  this  case  is  a  pure  supersonic  expansion,  the  exit  boundary  conditions  are  those  of  superson¬ 
ic  outflow;  i.e.,  the  exit  conditions  arc  determined  from  a  zeroth  order  extrapolation  of  the  dependent 
variables  at  they  - 1  grid  location.  Initial  conditions  for  this  case  were  the  uniform  distribution  of 
the  dependent  variables  as  determined  from  an  equilibrium  calculation  of  the  inflow  conditions. 

The  comparisons  of  the  NED  ANA  results  with  results  from  Molvik's  flow  solver  arc  pre¬ 
sented  in  Fig.  7  along  with  the  results  from  the  author’s  equilibrium  airflow  solver.  The  distribu¬ 
tions  of  the  species  and  the  flow  properties  are  nearly  indistinguishable  among  the  three  solutions, 
except  for  the  NO  comparison,  Fig.  7h.  Near  the  nozzle  exit,  where  the  density  has  dropped  by  an 
Older  of  magnitude  from  the  inflow  value,  the  NO  chemistry  appears  to  be  frozen. 
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c.  'Rxnperature  (normalized  by  T*) 
Figure?.  Super  duct. 
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d.  Internal  energy  (normalized  by  ei*) 


e.  Mixture  density  (normalized  by  p*) 


f.  Atomic  nitrogen  and  oxygen 
Figure  7.  Continued. 
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Figure  7.  Concluded. 
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4A  SUPERSONIC  DUCT  WITH  AREA  CHANGE  AND  NORMAL  SHOCK 

This  case  is  interesting  because  the  flow  downstream  of  the  normal  shock  is  subsonic. 
Consequently,  a  subsonic  boundary  condition  is  required  at  the  exit.  Initial  conditions  also  must 
reflect  this  fact.  Initial  conditions  were  determined  by  first  making  an  approximation  to  the  proper- 
.  ties  in  the  duct  from  a  perfect  gas  calculation.  Then  a  crude  equilibrium  air  calculation  was  made 
of  the  species  distribution  over  the  length  of  the  duct.  The  initial  shock  location  was  close  to  the 
expected  result.  All  of  the  dependent  variables  at  the  exit  were  determined  by  extrapolation  as 
explained  in  the  previous  subsection;  except  that  the  specific  total  energy  was  calculated  from  the 
exit  pressure  which  was  imposed  =  14.8  atm).  That  is, 


Knowing  the  mixture  temperature  and  the  constituents  at  the  exit,  NEQPAK  returns  the  specific  in- 
temal  energy.  The  total  specific  energy  is  then  found  by  adding  the  kinetic  energy,  as 

texit  =  f.7  +  JW''- 

The  comparisons  are  shown  in  Fig.  8.  Only  the  equilibrium  results  are  available  to  compare 
with  NEDANA,  because  the  version  of  Molvik’s  code  used  in  this  study  does  not  make  provisions 
for  subsonic  outflow  boundary  conditions.The  results  from  NEDANA  are  in  excellent  agreement 
with  the  equilibrium  air  solution  except  for  the  distribution  of  NO+  downstream  of  the  shock.  This 
results  because  the  equilibrium  air  solution  considers  eleven  species,  whereas  the  NEDANA 
computation  was  performed  with  a  seven  species  model. 

Note  the  changes  in  the  constituents  of  the  mixture  across  this  shock,  relative  to  the  changes 
discussed  with  regard  to  the  shock  and  contact  surface  in  the  shock  tube.  The  changes  in  pressure 
and  temperature  are  not  as  severe  in  the  duct  so  that  the  adjustment  process  is  not  as  smeared  as  in 
the  shock  tube.  The  constituents  change  cleanly  across  the  duct  shock,  except  for  a  spike  in  atomic 
oxygen.  It  is  also  important  to  recognize  that  in  the  duct,  substantial  dissociation  exists  on  both  sides 
of  the  shock;  whereas  in  the  shock  tube,  the  air  molecules  are  shocked  from  a  nearly  quiescent  state. 
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a.  Velocity  (normalized  by  a*) 


b.  Mixture  pressure  (normalized  by  p*) 


c.  Temperature  (normalized  by  T*) 
Figures.  Shock  tube  comparisons. 
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d.  Mixture  density  (normalized  by  ej*) 


e.  Mixture  density  (normalized  by  p*) 


f.  Atomic  nitrogen  and  oxygen 
Figures.  Continued. 


44 


Mole  Fraction  Mole  Fraction  Mole  Fraction 


AEDC-Tn-g4-18 


0.0008 

0.0006 

0.0004 

0.0002 

0 


2  4  6  8  10 

X 

h.  Nitric  oxide 


2  4  6  8  10 
X 


i.  Electrons 
Figure  8.  Concluded. 
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5.0  RESULTS  OF  TWO-  AND  THREE-DIMENSIONAL  COMPUTATIONS 

Extensive  flow  computations  were  performed  to  validate  and  demonstrate  the  methodology 
and  accuracy  of  the  two-dimensional  and  three-dimensional  modes  of  NEDANA.  The  classes  of 
problems  studied  include: 

1 .  Laminar  Separated  Flows  on  Hypersonic  Flat  PlateAVedges 

2.  Nonequilibrium  Flows  around  Hemisphere  Cylinders 

3.  Chimera  Domain  Decomposition  in  Nonequilibrium  Flows 

The  numerical  results  are  compared  to  both  experimental  data  and  numerical  results  from  existing 
flow  solvers  where  available. 

5.1  BOUNDARY  CONDITIONS 

In  this  section,  the  boundary  conditions  imposed  for  the  three  test  problems  are  discussed.  In 
all  cases,  the  flow  was  supersonic.  Therefore,  the  upstream  boundary  was  specified  with  free- 
stream  values.  The  flow  at  the  outflow  boundary  was  primarily  supersonic.  Only  a  small  subsonic 
region  existed  in  the  boundary  layer  at  the  outflow  boundary.  Therefore,  a  zeroth-order 
extrapolation  was  used.  Fluxes  at  symmetry  planes  were  calculated  by  reflecting  grid  points  across 
the  symmetry  boundary.  At  the  body  surface,  a  no-slip  condition  was  applied  such  that  all 
components  of  the  velocity  were  zero.  For  non-catalytic  walls,  a  zero  normal  gradient  in  the  species 
mass  fractions  was  set  at  the  wall.  For  fully  catalytic  walls,  the  species  mass  fractions  at  the  wall 
were  set  to  their  equilibrium  values  based  on  the  wall  temperature.  The  nonequilibrium  flow 
experiments  investigated  in  this  study  reported  a  constant  wall  temperature  of  555.55  K. 
Calculations  using  the  method  of  Prabhu  and  Erickson  (Ref.  48)  showed  that  for  this  temperature 
the  equilibrium  composition  of  air  was  equal  to  the  free-stream  composition.  Furthermore,  this  re¬ 
sult  was  nearly  independent  of  the  equilibrium  density.  Therefore,  the  species  mass  fractions  at  the 
wall  were  set  to  their  free-stream  values.  The  pressure  at  the  wall  was  set  assuming  zero  normal 
pressure  gradient.  The  temperatures  at  the  wall  were  set  to  the  reported  wall  temperature.  Park  (Ref. 
52)  reports  that  the  vibrational  and  electronic  temperatures  of  molecules  leaving  a  surface  are  nearly 
equilibrated  with  the  translational  wall  temperature.  Therefore,  the  translational-rotational  and 
vibrational-electronic  temperatures  at  the  wall  were  set  to  a  common  value. 
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Table  1.  Flow  Conditions  and  Geometries  for  Flat-PlateAVedge  Calculations 


Wedge  Angle  (deg) 

Re„(l/m) 

TJK) 

P^(Pa) 

TvallC^ 

15.05 

15.67 

4.885  X  (10)’’ 

40.31 

8.27 

294.4 

18.00 

15.58 

4.521  x(10)5 

42.53 

8.34 

294.4 

X,  m 


Figure  9.  Flat-plate/wcdgc  grid,  101  x  3  x  101. 

5.2  HYPERSONIC  LAMINAR  FLOW  OVER  A  FLAT-PLATE/WEDGE 

The  flow  over  a  compression  comer  fomied  by  the  intersection  of  a  flat  plate  and  a  wedge  was 
chosen  as  a  test  case  for  the  two-dimensional,  viscous  coding  of  NEDANA.  Measurements  for 
laminar,  attached  and  separated  flows  on  flat-plate/wedge  configurations  reported  in  the  CUBDAT 
database  by  Holden  (Ref.  53)  have  been  used  by  various  researchers  for  code  validation  (Refs.  54, 
55,  and  56).  Although  the  configuration  is  geometrically  simple,  the  physics  of  the  flow  field  are 
very  complicated  and  serve  as  an  ideal  test  case  for  viscous  hypersonic  phenomenon.  A  very  strong 
leading-edge  shock  is  generated  that  extends  downstream  and  intersects  the  wedge  flow  field.  In 
addition,  for  sufficiently  large  wedge  angles,  a  boundary-layer  separation  region  forms  in  the 
comer.  On  the  wedge,  downstream  of  this  separation  region,  the  flow  is  compressed  and  the 
boundary  layer  thins,  resulting  in  large  increases  in  skin  friction  and  heat  transfer.  Furthermore,  the 
compression  waves  produced  by  the  corner  coalesce  into  a  shock  wave  that  intersects  with  the 
leading-edge  shock,  generating  an  expansion  fan  and  shear  layer. 

The  CUBDAT  database  includes  pressure,  heat-transfer,  and  skin-friction  data  for  nominally 
two-dimensional  (finite-.span)  1 5.05-deg,  1  S-deg,  and  24-deg  wedges.  The  work  of  Rudy  et  al.  (Ref 
55)  concludes  that  the  inclusion  of  three-dimensional  effects  is  necessary  for  accurate  comparisons 
to  experimental  data  for  the  24-deg  wedge  case  due  to  the  larger  boundaiy-layer  .separation  region. 
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Therefore,  two-dimensional  comparisons  are  made  only  to  the  IS.OS-de^  and  I  ^-deg  wedge  data. 
The  conditions  for  the  simulations  are  listed  in  Table  1.  Note  that  the  Reynolds  number  was  low 
enough  that  the  flow  remained  laminar.  In  addition,  the  high  Mach  number  is  achieved  by  an 
extremely  low  free-stream  temperature.  Therefore,  the  flow  is  assumed  a  perfect  gas  and  not  a 
chemically  reacting  gas.  As  a  result,  the  flat-plate  and  wedge  data  comparisons  allow  for  a 
validation  of  the  NEDANA  flow  solver  in  the  absence  of  the  additional  complexities  of  turbulence 
and  chemistry. 

For  both  wedge  angles,  the  length  of  the  flat  plate  and  wedge  are  0.4396  m  and  0,3048  m, 
respectively.  The  test  geometries  were  modeled  as  two-dimensional  compression  comers  with 
101  X  3  X  101  grids.  (Note  that  the  current  version  of  NEDANA  is  a  three-dimensional,  finite- 
volume  flow  solver  with  the  grid  system  described  in  Appendix  6.  Therefore,  a  minimum  of  three 
grid  points  is  required  in  all  coordinate  directions.)  A  viscous  wall  spacing  of  1.0  x  10'^  m  was 
used.  This  viscous  spacing  produced  values  of  <  0.1 .  The  calculation  of  y*  is  useful  for  deter¬ 
mining  the  minimum  grid  spacing  required  for  accurate  viscous  solutions  even  in  the  absence  of 
turbulence.  A  typical  grid  is  shown  in  Fig.  9.  The  current  version  of  NEDANA  has  been  devel¬ 
oped  for  nonequilibrium  computations  and,  therefore,  receives  all  thermodynamic  and  transport 
properties  from  NEQPAK.  As  a  consequence,  the  perfect-gas  flow  was  modeled  as  a  chemically 
frozen  mixture  of  Wj  and  O2,  representing  air.  This  fact  is  noteworthy  because  the  NEQPAK 
curve  fits  for  thermodynamic  and  transport  properties  were  formulated  for  100  K  to  30,000  K. 

Steady-state  computations  were  made  with  the  NEDANA  flow  solver  for  both  the  15.05-^eg 
and  1 6-deg  wedge  configurations  using  a  global  time  step  and  a  maximum  CFL  number  of  five.  The 
dissipation  parameters  ((see  Eq.  (S2))  were  set  at  Kj  =  0.8  and  iq  =  1 .8.  Numerical  tests  showed  that 
the  use  of  local  time  stepping  greatly  increased  the  number  of  iterations  for  convergence.  This  find¬ 
ing  may  be  due  to  the  unsteady  nature  of  the  separated  flow  region.  Steady-state  solutions  were 
typically  achieved  in  10,000  iterations  with  a  reduction  in  residual  of  four  orders  of  magnitude. 
Pressure,  heat-transfer,  and  skin-friction  results  for  the  15.05-<feg  wedge  case  are  presented  in  Figs. 
10, 1 1 ,  and  12.  The  results  are  compared  to  experimental  data.  Solutions  were  also  computed  using 
Molvik’s  TLIFF  (Ref.  49)  flow  solver  and  the  AEDC  flow  solver  XMERA  (Ref.  42).  These  results 
are  included  to  assess  the  performance  of  NEDANA  versus  other  state-of-the  ait  flow  solvers.  Re¬ 
sults  from  all  three  flow  solvers  are  in  good  agreement  with  the  experimental  data.  Note,  the  results 
of  all  three  flow  solvers  may  differ  slightly  due  to  the  use  of  different  curve  fits  for  thermodynamic 
and  transport  properties. 


48 


AEDC-TR-94-ie 


Figure  10.  Pressure  distributions  for  laminar  flow  over  a  15.05-deg  wedge, 
=  15.67. 


Figure  11.  Heat-transfer  distributions  for  laminar  flow  over  a  15.05*deg  wedge, 
=  15.67. 
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Figure  12.  Skin  friction  distributions  for  laminar  flow  over  a  lS.OS-deg  wedge, 
M„=  15.67. 


The  calculation  of  the  size  of  the  sq)aration  region  that  fonns  at  the  juncture  of  the  flat  plate/ 
wedge  is  extremely  sensitive  to  the  numerical  dissipation  of  the  flow  solver.  A  flow  solver  with  low 
numerical  dissipation  should  accurately  predict  both  the  size  and  shape  of  the  separation  region.  As 
the  dissipation  of  the  numerical  scheme  increases,  the  size  of  the  predicted  separation  region 
conversely  decreases.  Extremely  dissipative  schemes  may  actually  suppress  the  formation  of  a 
separation  region  entirely.  The  effect  of  the  separated  flow  is  to  decrease  both  the  skin  friction  and 
heat  transfer  at  the  wall.  As  a  result,  the  extent  of  the  separation  region  is  marked  by  a  drop  in  the 
skin-friction  and  heat-transfer  distributions.  With  this  phenomenon  in  mind,  an  examination  of  the 
skin-friction  and  heat-transfer  distributions  reveals  that  NEDANA  and  XMERA  are  slightly  more 
dissipative  than  TUFF.  This  result  was  anticipated  because  the  TUFF  flow  solver  is  based  on  an 
upwind  algorithm,  whereas  the  NEDANA  and  XMERA  flow  solvers  are  based  on  central  differ¬ 
ence  schemes  with  constant  numerical  dissipation  parameters. 

To  investigate  the  effects  of  decreased  numerical  dissipation,  a  series  of  computations 
were  performed  with  NEDANA  in  which  the  dissipation  parameters  were  lowered  to  =  0.4  and 
=  0.8.  Additional  computations  were  performed  where  the  dissipation  parameters  were  also 
scaled  by  the  normalized  velocity  such  that  IC2  — >  0.0  and  »  0.0  at  the  no-sIip  wall.  This 
scaling  was  performed  to  further  lower  the  dissipation  in  the  viscous  region  near  the  wall.  The 
decreased  numerical  dissipation  did  not  significantly  improve  the  NEDANA  results.  Only  a 
slight  increase  in  the  size  of  the  separation  region  was  achieved,  and  at  a  cost  of  decreased 
numerical  stability.  A  further  numerical  experiment  was  performed  where  the  viscous  spacing 
was  increased  to  1  .Ox  10*^  m.  The  ensuing  computation  completely  failed  to  predict  the  separation 
region.  This  result  would  indicate  that  adequate  viscous  grid  resolution  is  more  significant  than 
the  magnitude  of  the  dissipation  parameters  employed. 
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Similar  results  are  presented  in  Figs.  13, 14,  and  15  for  the  1  ^-deg  wedge  case.  The  NEDANA 
and  XMERA  flow  solvers  performed  similarly.  However,  the  upwind  Roe  scheme  of  TUFF  initially 
generated  a  nonphysical  solution  in  the  vicinity  of  the  separation  region  where  the  flow  switches 
from  supersonic  to  subsonic.  This  nonphysical  behavior  is  common  to  upwind  Roe  schemes  and 
can  be  eliminated  by  applying  an  entropy  fix  (eigenvalue  limiter)  to  the  eigenvalues  (Ref.  49).  By 
increasing  the  absolute  value  of  the  entropy  fix,  physically  consistent  solutions  were  obtained. 
However,  the  TUFF  results  then  offered  no  advantage  over  the  central  difference  NEDANA  and 
XMERA  flow  solvers.  AH  three  flow  solvers  failed  to  completely  predict  the  extent  of  the  separa¬ 
tion  region. 


Figure  13.  Pressure  distributions  for  iaminar  flow  over  an  18-deg  wedge, 
M„  =  15.58. 


The  results  of  the  hypersonic  laminar  flat-plate/wedge  cases  have  shown  NEDANA  to 
perform  as  well  as  current  state-of-the  art  flow  solvers.  The  NEDANA  flow  solver  proved  to  be 
more  stable  and  to  converge  more  rapidly  for  h3rpersonic  flow  conditions  than  the  XMERA  flow 
solver.  Where  the  NEDANA  flow  solver  was  able  to  run  with  a  CFL  number  of  five  and  converge 
in  10,000  iterations,  the  XMERA  flow  solver  was  limited  to  a  CFL  number  of  one-tenth  and 
required  ](X),000  iterations  to  converge.  Where  the  NEDANA  code  proved  to  be  robust  for  both 
wedge  geometries,  the  TUFF  flow  solver  was  susceptible  to  spurious  solutions,  even  for  small 
changes  in  wedge  angle. 
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Figure  14.  Heat-transfer  distributions  for  laminar  flow  over  an  18*deg  wedge, 
=  15.58. 


Figure  15.  SIdn-fIriction  distributions  for  laminar  flow  over  an  18-deg  wedge, 
M„=  15.58. 
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5.3  NONEQUILIBRIUM  FLOW  AROUND  HEMISPHERE  CYLINDERS 

The  computation  of  the  thermal  and  chemical  nonequilibrium  flow  fields  around  hemisphere 
cylinders  (Refs.  57  and  58)  was  selected  as  a  validation  case  for  the  nonequilibrium  capabilities  of 
NEDANA.  Computations  were  made  for  the  three  geometries  listed  in  Table  2.  The  three 
hemisphere  cylinders  have  different  radii  and  are  referred  to  as  model  1,  model  2,  and  model  3.  The 
nominal  flow  conditions  for  all  three  models  arc  =  9.8,  Re^  -  3.025  x  (10)^  (I/m),  T_^  «  450  K, 
p  =  230  Pa,  and  T^,  =  555.55  K.  The  conditions  are  such  that  the  flow  is  in  chemical  and  thermal 
nonequilibrium.  The  frce-stream  total  enthalpy,  is  9.14  MJfkg.  A  five  species  air  chemistty, 
two-temperature  model  was  applied  using  the  reaction  rates  of  Park  (Ref.  36).  In  this  work,  the 
dissociation  reactions  arc  governed  by  a  generic  temperature  where  a  =  0.3  (see  Eq.  (31)).  Heat- 
transfer  measurements  were  obtained  with  both  thin-film  and  thermocouple  gauges.  The  thin-film 
gauges  were  coated  with  mother  of  pearl  to  minimize  wall  catalycity,  whereas  the  thermocouple 
gauges  allowed  for  wall  catalysis.  Models  1  and  2  were  fitted  with  stagnation  point  thin-film  gauges 
for  measuring  heat  transfer,  and  model  3  was  fitted  with  a  stagnation  point  thermocouple  gauge  for 
measuring  heat  transfer. 

Table  2.  Hemisphere/Cylinder  Geometries  and  Stagnation  Point 
Catalytic  Boundary  Conditions 


Model 

Rn  (m) 

L(ffi) 

Wall  B.C. 

1 

0.0127 

0.1143 

Noncatalytic 

2 

0.0254 

0.2032 

Noncatalytic 

3 

0.0508 

0.4064 

Catalytic 

Grid  resolution  studies  were  performed  with  the  current  flow  solver  on  the  intermediate¬ 
sized  model,  model  2,  to  determine  the  grid  requirements  for  achieving  valid  pressure  and  heat- 
transfer  results  for  nonequilibrium  flows.  Josyula  and  Shang  (Ref.  57)  reported  using  a  mesh 
system  with  40  points  in  the  body-tangential  direction  and  50  points  in  the  body-normal  direction 
with  a  viscous  wall  spacing  of  5.0  x  (10)'^  m.  Using  the  results  of  Josyula  and  Shang  as  a  guide, 
the  system  of  meshes  listed  in  Table  3  was  developed.  A  typical  grid  is  shown  in  Fig.  16.  Grids  I, 
2,  and  3  were  developed  such  that  the  grid  distribution  functions  in  the  body-tangential  and  body- 
normal  directions  were  the  same.  This  criterion  was  enforced  to  maintain  the  same  relative 
resolutions  in  all  three  grids.  Grid  4  was  created  with  the  same  number  of  grid  points  as  Grid  1  but 
with  the  viscous  spacing  of  Grid  3.  Grid  4  was  developed  to  assess  the  relative  importance  of  the 
viscous  spacing  and  the  total  number  of  points  in  the  body-normal  direction.  (Grids  5, 6,  and  7  will 
be  discussed  later.)  Note  that  NEDANA  was  used  in  a  fully  three-dimensional,  finite-volume 
mode  with  the  grid  system  described  in  Appendix  B.  Therefore,  five  grid  points  were  required  in 
the  body-azimuthal  direction. 
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Table  3.  Details  of  Grid  Resolution  Study 


Grid 

Dimensions 

No.  Points 

An(ffi) 

Iteration  for  Conversion 

1 

51  x5x51 

13005 

1.0  X  (10)-^ 

1100 

2 

75  X  5  X  75 

28125 

7,0  X  (10)-^ 

1900 

3 

101  x5x  101 

51005 

5.0  X  (10)-* 

2850 

4 

51  x5x51 

13005 

5.0  X  (10)-^ 

1300 

5 

51  X  5  X  51 

13005 

1.0  X  (10)-^ 

3700 

6 

51  X  5x75 

19125 

1.0  X  (10)-^ 

4500 

7 

51  x5x75 

19125 

5.0  X  (10)  '' 

5000 

Steady-state  computations  were  made  with  NEDANA  for  Grids  1,  2,  3,  and  4  using  local  time 
stepping  and  a  maximum  CFL  number  of  ten.  The  dissipation  parameters  were  set  at  =  0.8  and 
=  1.8  for  all  cases.  Numerical  te.sts  showed  that  the  use  of  a  local  time  step  gave  the  optimal 
convergence  rate  for  the  hemisphere  cylinder  flow  fields.  (This  result  is  contrary  to  the  findings  of 
the  flat-plate/wedge  investigation  and  emphasizes  the  uncertainties  of  using  local  time-stepping 
methods.)  The  solutions  were  considered  converged  when  the  heat  transfer  to  the  wall  varied  by 
less  than  1  percent  over  a  characteristic  time  step  based  on  the  minimum  time  step  for  the  grid. 
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where  =Rn/  u^.  The  number  of  time  steps  required  to  obtain  a  steady-state  solution  is  also  listed 
for  each  grid  in  Table  3.  The  NEDANA  flow  solver  operated  at  a  speed  of  2.03  x  10*^  cpu-s^t/iter 
on  a  Convex  3840  computer  for  five  species  in  chemical  and  thermal  nonequilibrium. 

Surface  pressure  and  heat-transfer  distributions  for  the  grid  resolution  study  of  model  2  with 
grids  1,  2,  3,  and  4  are  presented  in  Figs.  17  and  18.  The  comparisons  with  experimental  measure¬ 
ments  for  pressure  are  excellent.  Furthermore,  the  calculation  of  surface  pressure  is  largely 
independent  of  the  grid  used.  The  pressure  distributions  for  all  grids  differed  by  less  than  2  percent. 
The  comparisons  with  experimental  heat-transfer  measurements  show  a  greater  grid  dependency. 
As  the  grid  resolution  and  viscous  spacing  are  refined,  the  computed  heat  transfer  increases,  llie 
viscous  spacing  normal  to  the  wall  is  critical  in  predicting  surface  heat  transfer.  This  fact  is  evident 
in  the  comparison  of  the  surface  heat  transfer  for  grids  3  and  4.  Grid  4  has  the  same  initial  viscous 
spacing  as  grid  3,  but  only  half  the  number  of  points  normal  to  the  body.  Nevertheless,  the  surface 
heat-transfer  distribution  obtained  on  grid  4  is  in  better  agreement  with  the  results  obtained  on  grid 
3  than  those  obtained  on  grids  1  or  2.  A  comparison  of  the  stagnation-point  heat  transfer  for  grids 
1,  2,  and  4,  shows  they  differ  from  that  of  grid  3  by  approximately  20  percent,  10  percent,  and  5 
percent,  respectively.  The  comparisons  are  taken  relative  to  grid  3,  because  it  is  assumed  the  grid 
with  the  best  resolution  will  produce  the  most  accurate  results.  Clearly,  the  computed  surface  heat- 
transfer  distributions  are  grid  dependent.  For  grid  3,  y'''  =  0.3  at  the  stagnation  point,  and  down¬ 
stream  of  the  spherical  surface  y*  <^\  .5.  A  viscous  resolution  corresponding  to  <  1  was  initially 
expected  to  yield  grid-independent  results.  However,  the  criteria  for  the  accurate  calculation  of 
surface  heat  transfer  proved  to  be  more  complex. 

Siddiqui  et  al.  (Ref.  59)  report  similar  difficulties  in  computing  grid-independent  surface  heat 
transfer  using  a  variety  of  thin-layer  Navier-Stokes  algorithms.  Siddiqui  et  al.  computed  the  flow 
over  spherically  blunt  cones  at  similar  Mach  numbers  and  Reynolds  numbers.  (However,  the  free- 
stream  temperature  and  velocity  were  such  that  the  flows  were  not  chemically  reacting.)  They  de¬ 
veloped  an  expression  to  estimate  the  required  viscous  spacing  as  a  function  of  Mach  number  and 
Reynolds  number.  Applying  this  expression  to  the  hemisphere  cylinders  of  this  study,  the  estimated 
viscous  spacing  required  for  grid-independent  heat  transfer  is  5.0  x  (I0‘^)  m.  To  generate  agrid  with 
this  viscous  spacing  and  maintain  a  similar  grid  distribution  as  grids  1, 2  and  3  is  not  practical.  The 
resulting  grid  would  require  hundreds  of  points  normal  to  the  body.  Also,  note  that  the  number  of 
iterations  to  convergence,  as  well  as  the  cpu-s/iter  would  increase  beyond  practicality.  Hierefore,  a 
compromise  was  made  by  holding  the  number  of  points  normal  to  the  body  constant.  As  a  result,  a 
decrease  in  the  viscous  spacing  produces  poorer  resolution  of  the  shock  layer.  This  reduced  resolu¬ 
tion  in  the  shock  layer  is  of  importance  in  nonequilibrium  flows  where  the  chemical  reactions  can 
affect  the  heat  transfer  to  the  surface.  This  need  to  accurately  resolve  the  chemical  reactions  in  the 
shock  layer  is  an  additional  constraint  not  imposed  on  the  perfect-gas  work  of  Siddiqui  et  al. 
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To  assess  the  effects  of  even  smaller  viscous  spacings,  a  second  series  of  grid  resolution 
computations  was  performed  with  viscous  spacings  approaching  those  recommended  by  Siddiqui 
et  al.  These  grids  are  listed  in  Table  3  as  grids  5, 6,  and  7.  Surface  pressure  and  heat-transfer  distri¬ 
butions  for  grids  5,  6,  and  7  are  presented  in  Figs.  19  and  20.  The  comparisons  with  experimental 
measurements  for  pressure  are  once  again  excellent.  Furthermore,  the  calculation  of  surface 
pressure  is  largely  independent  of  the  grid  used.  The  pressure  distributions  for  all  grids  differed  by 
less  than  2  percent.  The  pressure  distributions,  as  expected,  also  are  in  agreement  with  the  results 
obtained  from  grids  1-4.  The  comparisons  with  experimental  heat-transfer  measurements  for  grids 
5,  6,  and  7  are  similar.  However,  the  heat-transfer  distributions  still  exhibit  a  sensitivity  to  the 
viscous  spacing.  The  stagnation-point  values  of  heat  transfer  for  grids  S  and  6  differ  by 
approximately  20  and  10  percent,  respectively,  from  those  of  grid  7.  For  grid  7,  «  0.10  at  the 

stagnation  point.  Downstream  of  the  spherical  surface,  =  0.15.  For  grids  5  and  6,  “  0.15  at  the 

stagnation  point.  Downstream  of  the  spherical  surface,  y*  «  0.20.  These  results  further  demonstrate 
the  difficulties  in  achieving  truly  grid-independent,  heat-transfer  predictions.  However,  these 
results  coupled  with  the  results  of  Sec.  5.2  also  demonstrate  that  for  the  cases  studied,  grids  with 
viscous  spacings  corresponding  to  y*  »  0.1  produced  heat-transfer  predictions  within  the 
experimental  scatter. 


Figure  19.  Pressure  distributions  for  hemisphere/cylinder  model  2  for  grids  5, 6,  and  7. 
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Figure  20.  Heat-transfer  distributions  for  hemisphere/cylinder  model  2  for 
grids  5,  6,  and  7. 

The  effects  of  viscous  resolution  on  heat-transfer  predictions  have  been  well  documented. 
However,  the  importance  of  resolving  the  shock  layer  should  not  be  overlooked.  This  constraint  is 
seen  in  the  comparisons  of  surface  heat  transfer  for  grids  5  and  6.  Both  grids  have  the  same  initial 
viscous  spacing.  They  only  differ  in  the  number  of  points  and  distribution  of  points  normal  to  the 
body.  The  computed  stagnation-point  heat  transfer  for  these  two  grids  differs  by  approximately  8 
percent,  demonstrating  the  importance  in  resolving  the  shock  layer.  Also,  the  heat-transfer 
distributions  for  grids  S,  6,  and  7  differ  from  those  of  grids  1-4  downstream  of  the  hemisphere  nose. 
However,  instead  of  predicting  an  increased  heat-transfer  rate  to  the  surface,  they  show  a  lower 
heat-transfer  rate  to  the  body.  This  trend  is  more  in  line  with  the  experimental  data.  Siddiqui  et  al. 
(Ref.  59)  state  that  this  lower  heat-transfer  rate  is  due  to  a  better  resolution  of  the  expansion  waves 
emanating  from  the  spherical  region.  This  consequence  of  better  resolving  the  expansion  region  of 
the  flow  further  emphasizes  the  difficulties  associated  with  accurate  prediction  of  heat-transfer  data. 

Based  on  the  results  of  the  two  grid  resolution  studies,  the  geometries  of  models  1,  2,  and  3 
were  modeled  with  51  x  5  x  75  grids  with  initial  viscous  spacings  of  5.0  x  (lO)'^m.  These  viscous 
spacings  produced  values  of  y** »  0. 1 .  Surface  pressure  and  heat-transfer  distributions  for  models  1 , 
2,  and  3  are  presented  in  Figs.  21-26.  The  surface  pressure  distributions  are  in  excellent  agreement 
with  the  experimental  data  for  all  three  models.  The  agreement  with  the  stagnation-point  heat- 
transfer  data  is  reasonable  for  all  three  models.  The  agreement  with  the  heat-transfer  data  on  the 
afterbodies  is  within  the  experimental  data  scatter.  For  model  3,  both  a  fully  catalytic  wall  and  a 
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noncatalytic  wall  boundary  condition  were  applied.  The  unc^ainty  in  the  stagnation-point  heat- 
transfer  measurement  of  model  3  is  reported  as  zero  because  only  one  measurement  was  taken.  Note 
that  the  fiilly  catalytic  wall  increases  the  stagnation-point  heat  transfer  by  58  percent  over  that  of 
the  noncatalytic  wall.  Generally,  the  model  surface  will  exhibit  finite  catalytic  rates.  Therefore,  the 
true  heat-transfer  distribution  is  bounded  by  the  fully  catalytic  and  noncatalytic  values.  Clearly,  the 
accurate  modeling  of  wall  catalycity  represents  a  significant  challenge  in  the  accurate  prediction  of 
surface  heat  transfer  in  both  computational  and  experimental  simulations. 


Figure  21.  Pressure  distributions  for  hemisphere/cylinder  model  1. 
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Figure  22.  Heat*transfer  distributions  for  hemisphere/cylinder  model  1. 
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Figure  26.  Heat-transfer  distributions  for  taemisphere/cylinder  model  3. 

Temperature  distributions  along  the  stagnation  streamline  for  models  1,  2,  and  3  are 
presented  in  Fig.  27.  The  conditions  directly  behind  the  shock  are  determined  by  the  free-stream 
conditions.  Therefore,  the  peak  post-shock  temperatures  for  all  three  models  are  nearly  the  same. 
(Small  differences  in  the  post-shock  temperatures  may  be  attributed  to  differences  in  grid  resolu¬ 
tion.)  Note,  the  translational  temperature  rises  rapidly  behind  the  shock.  However,  a  much  longer 
time  is  required  for  the  vibrational  modes  of  the  molecules  to  become  excited.  As  a  result,  the 
vibrational  temperature  lags  the  translational  temperature,  producing  a  state  of  thermal  nonequilib¬ 
rium.  The  time  or  distance  that  is  required  for  the  translational  and  vibrational  temperatures  to 
equilibrate  is  dependent  on  this  initial  departure  from  equilibrium,  and  is  therefore  similar  for  all 
three  models.  However,  the  shock-stand  off  distance  for  each  model  is  proportional  to  the  nose 
radius.  Therefore,  as  the  nose  radius  decreases,  a  larger  extent  of  the  shock  layer  is  in  a  state  of 
thermal  nonequilibrium.  Also,  note  that  the  peak  post-shock  temperature  occurs  closer  to  the  body 
with  decreasing  nose  radius.  This  condition  produces  steeper  temperature  gradients  near  the  body, 
leading  to  increased  heat-transfer  rates.  These  findings  are  consistent  with  hypersonic  boundary- 
layer  theory  (Ref.  60)  which  states  that  q<K\ljRn . 
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Figure  27.  Temperature  distributions  along  stagnation  streamline  for 
models  1, 2,  and  3. 


5.4  CHIMERA  DOMAIN  DECOMPOSITION 


As  a  demonstration  of  the  chimera  capability  of  NEDANA,  a  computation  for  the  hemisphere 
cylinder  model  2  was  made  with  the  computational  domain  defined  by  two  overlapping  grids 
consisting  of  45  x  5  x  75  points  on  the  forebody  and  25  x  5  x  5 1  points  on  the  afterbody.  The  two- 
grid  system  is  shown  in  Fig.  28.  The  two  grids  had  different  numbers  of  points  and  grid  distributions 
in  both  directions.  Therefore,  there  was  no  exact  overlap  (direct  injection)  of  points  between  the  two 
grids.  The  one  grid  constraint  that  was  applied  was  that  the  first  point  off  the  wall  for  the  two  grids 
has  the  same  viscous  spacing.  This  constraint  was  applied  because  of  the  extreme  dependence  of 
heat-transfer  results  on  viscous  .spacing.  To  facilitate  the  comparisons  of  the  results  obtained  on  a 
single  grid,  the  two  chimera  grids  were  given  initial  viscous  spacings  equal  to  that  of  grid  1  in  Table 
3.  Hence,  comparisons  are  also  made  to  the  re.sults  of  grid  1.  Surface  pressure  and  heat-transfer 
distributions  for  grid  1  and  the  chimera  grid  system  are  compared  in  Figs.  29  and  30.  The 
distributions  for  the  chimera  grids  are  almost  identical  to  tho.se  of  grid  1  even  with  the  mismatches 
in  grid  resolution.  The  pressure  and  heat-Uansfer  distributions  are  also  smooth  in  the  overlap  region 
between  the  forebody  and  afterbody  grids. 
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Figure  28.  Chimera  hemisphere/cylinder  grid,  45  x  5  x  75  and  25  x  5  x  51. 
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Figure  29.  Pressure  distributions  for  model  2  with  single  and  chimera  grid  system. 
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Figure  30.  Heat-transfer  distributions  for  model  2  with  single  and 
chimera  grid  system. 

These  results  demonstrate  the  ability  of  the  chimera  grid  system  to  reproduce  the  surface  pres¬ 
sure  and  heat-transfer  results  of  the  single  grid.  However,  questions  have  been  raised  about  the  ap¬ 
plicability  of  the  chimera  scheme  to  hypersonic  nonequilibrium  flows,  in  particular,  the  ability  of 
the  scheme  to  capture  strong  shocks,  as  well  as  conserve  species  concentrations  in  the  vicinity  of 
grid  overlaps.  In  Figs.  31  —  34,  contour  plots  are  presented  for  the  translational-rotational  and 
vibrational-electronic  temperatures,  and  for  the  mass  fractions  of  diatomic  and  monatomic  oxygen. 
Contour  lines  are  plotted  on  both  grids  in  the  overlap  region.  Therefore,  any  discrepancies  in  the 
temperatures  or  mass  fractions  between  the  grids  should  appear  as  mismatches  in  the  contour  lines. 
Small  discrepancies  are  perceivable  in  the  overlap  region.  They  occur  primarily  in  the  plots  of 
translational-rotational  temperature  and  not  in  the  plots  of  vibrational-electronic  temperature  or 
mass  fractions. 
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Figure  31.  Translational-rotational  temperature  contours  for  chimera  grid  system. 


Figure  32.  Vibrational-electronic  temperature  contours  for  chimera  grid  system. 
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The  contour  plots  provide  a  means  for  examining  the  qualitative  properties  of  the  chimera 
scheme  in  the  overlap  region.  However,  contour  plots  may  be  misleading  because  they  depend  on 
the  contour  levels  chosen.  Therefore,  temperature  and  mass  fraction  profiles  normal  to  the  body  at 
the  location  indicated  by  the  arrows  (x  »  0.036)  are  presented  in  Figs.  35  and  36.  The  profiles 
extracted  tom  the  fotebody  and  afterbody  grids  are  nearly  identical.  Also,  note  the  effects  of  both 
thermal  and  chemical  nonequilibrium.  Ifie  overlap  region  and  the  profile  location  extend  tom  the 
free  stream  through  the  shock  layer  to  the  surface  of  the  cylinder.  Examining  Fig.  35,  the  transla¬ 
tional-rotational  and  vibrational-electronic  temperatures  are  seen  to  be  in  equilibrium  in  the  toe 
stream.  As  the  distance  to  the  body  decreases,  the  profiles  correspond  to  a  fluid  element  that  crossed 
the  oblique  shock  just  upstream  of  the  profile  station.  Therefore,  the  translation-rotational  temper¬ 
ature  increases  while  the  vibrational-electronic  temperature  lags.  This  behavior  is  similar  to  the 
behavior  along  the  stagnation  line  that  was  discussed  in  the  previous  section.  The  profiles  at  the 
stagnation  line,  by  definition,  correspond  to  a  single  streamline  and  shock  crossing  position. 
However,  unlike  the  profile  along  the  stagnation  line,  the  current  profiles  correlate  to  many  stream¬ 
lines  and  shock  crossing  positions.  As  the  distance  to  the  body  in  Fig.  35  decreases  fiirther,  the  pro¬ 
file  corresponds  to  fluid  elements  that  have  come  through  the  much  stronger  shock  of  the  stagnation 
region.  These  fluid  elements  have  also  been  rapidly  expanded  around  the  hemisphere  nose  and  are, 
to  a  great  extent,  vibrationally  and  chemically  frozen.  An  examination  of  the  temperature  profiles 
near  the  surface  show  that  the  vibrational-electronic  temperature  is  frozen  at  a  much  higher  level 
than  the  translational-rotational  temperature. 

The  effects  of  chemical  nonequilibrium  and  freezing  are  seen  in  the  mass  fraction  profiles  of 
Fig.  36.  The  increase  in  the  translational-rotational  temperature  that  is  a  result  of  crossing  the 
oblique  shock  is  not  accompanied  by  a  change  in  the  mass  fractions  of  oxygen  as  is  experienced  at 
the  stagnation  line.  This  behavior  is  a  result  of  the  weaker  oblique  shock.  The  oblique  shock 
produces  lower  temperature,  and  density  rises  relative  to  the  shock  at  the  stagnation  region.  The 
dissociation  rates  behind  the  oblique  shock  are,  therefore,  much  slower.  However,  significant 
variations  in  the  mass  fractions  are  seen  as  the  distance  to  the  body  decreases  further.  Near  the  body, 
the  profiles  correspond  to  fluid  elements  that  have  passed  through  the  much  stronger  shock  of  the 
stagnation  region  and  the  subsequent  expansion  region.  Therefore,  the  composition  of  the  gas  is 
partially  frozen.  This  effect  of  chemical  nonequilibrium  is  evident  in  the  increased  degree  of 
dissociation  of  diatomic  oxygen  near  the  body.  If  the  gas  were  in  equilibrium  at  the  local 
translational  temperature,  the  degree  of  dissociation  of  the  gas  would  be  decreasing  with  decreasing 
temperature  near  the  body. 
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Figure  35.  Temperature  proflles  in  chimera  overlap  region. 


Figure  36.  Mass  fraction  Oj  and  O  profiles  in  chimera  overlap  region. 
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6.0  CONCLUSION 

A  three-dimensional  nonequilibrium  chimera  flow  solver  was  developed  and  tested  against 
existing  flow  solvers  and  experimental  data.  The  present  flow  solver,  NEDANA,  uses  an  implicit, 
time-marching  algorithm  that  is  central  di^erenced  and  made  TVD  (Total  Variation  Diminish¬ 
ing)  through  flux  limiters.  The  LIM  (Locally  Implicit  Method)  scheme  is  used  at  each  time  step 
to  obtain  a  solution  to  the  nonlinear  set  of  equations.  Flexibility  in  the  use  of  various 
nonequilibrium  chemistry  models  is  made  possible  by  the  availability  of  an  AEDC-developed 
nonequilibrium  chemistry  package,  NEQPAK.  Use  of  such  a  chemistry  package  makes 
NEDANA  unique  in  the  field. 

One-dimensional  comparisons  with  the  results  from  state-of-the  ait  flow  solvers  woo 
excellent.  It  was  expected  at  the  outset  that  sharper  shocks  in  one-dimensional  flows  would  be 
traded  for  a  simpler  algorithm,  but  this  was  not  the  case.  The  simpler  algorithm  (NEDANA) 
captured  shocks  as  well  as  those  of  the  more  complex  solvers.  However,  it  is  also  evident  that 
NEDANA  does  not  handle  contact  surfaces  as  well.  But,  these  surfaces  only  occur  in  unsteady 
flows,  which  are  of  somewhat  less  interest  at  AEDC.  Hie  three-dimensional  NEDANA  flow  solver 
has  been  extensively  tested  against  other  flow  solvers  and  available  experimental  data  for  perfect- 
gas  and  thermo-chemical  nonequilibrium  air  chemistry.  As  expected,  the  NEDANA  code  proved  to 
be  slightly  more  dissipative  than  upwind  schemes,  but  the  effects  of  increased  dissipation  were 
marginal.  A  nonequilibrium  chimera  capability  was  also  demonstrated,  emphasizing  the  flexibility 
of  the  NEDANA  methodology. 

The  present  solver’s  simplicity  helps  achieve  the  goal  of  having  solvers  that  are  relatively  easy 
to  use,  understand,  and  modify. 
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The  purpose  of  this  appendix  is  to  begin  with  very  few  simplifying  assumptions,  and  to 
develop  a  rather  general  set  of  conservation  equations  modeling  viscous  nonequilibrium  flow.  Then, 
various  assumptions  are  made  to  obtain  a  simplified,  reduced  equation  set  from  the  more  general. 


First,  the  internal  energy  per  unit  mass,  of  a  given  species,  s,  is  assumed  to  be  composed 

as  follows: 


Cj  —  "I"  ■I'  "I" 


s  =  1,..  .,n^ 


where  e,,,  e^jUnd  are  translational,  rotational,  vibrational,  and  electronic  parts,  respectively. 
Also,  ru  is  the  total  number  of  species.  Each  of  the  above  energies  has  an  associated  temperature. 
Specifically,  in  thermodynamic  equilibrium,  the  most  probable  fractions,  q  =  t.r,  v,  e.  of  the 
particles  of  a  given  species,  i,  at  the  discrete  energies,  ,  q  =  t,r,v,  e,  respectively,  are  described 
by  Boltzmann  distributions, 


,s  =  i,  ...,ns 


with  temperatures,  T^,  q  =  t,r,  v,  e,  respectively.  Here,  i  is  an  index  over  discrete  energy  states. 
Also,  d^l  is  the  number  of  distinct  slates  with  energy,  ,  and  is  the  partition  function. 


J 


The  Boltzmann  constant  is  k  ='R/Na,  where  7^is  the  universal  gas  constant  and  is  Avagadro’s 
number.  To  obtain  the  four  energies  for  the  mixture,  a  mass  average  is  performed: 


^7 


q  =  t’,e. 


Here,  is  the  partial  density  of  species,  s,  and  p  is  the  mixture  density  given  by  the  sum  of  the 
partial  densities.  Also,  the  mixture  internal  energy  is  given  by: 


77 


AEDC>TR-94-ie 


In  contrast  to  the  species  that  have  an  internal  structure,  the  free  electron  has  only  one  form  of 
energy,  namely  translational  energy.  Specifically, 

fe-  =  «r,c-  =  =  0.  (A-1) 

Also,  because  elemental  species  do  not  have  components  that  vibrate  with  respect  to  each  other, 

«i'..  =  0,  molecule.  (A-2) 


Similarly,  the  enthalpy  per  unit  mass,  h„  for  a  given  species,  s,  is  assumed  to  be  the  sum  of 
translational,  rotational,  vibrational,  and  electronic  parts; 

ft,  =  ht,,  +  hr,,  +  hv.,  +  fte.».  »=  1 . 


These  are  related  to  the  internal  energies  according  to 


ft 


9,*  ~ 


I 


+ 


'RTt,, 

M, 


q  =  i 


(A-3) 


where. Vi,  is  the  molecular  weight  of  species,  s.  Thus,  by  Eq.  (A-1),  the  free  electron  has  only  one 
form  of  enthalpy: 

ft,-  —  h.fg-„  hj.g~  =  ft„_,—  =  hf^f-  —  0.  (A-4) 


Also,  by  Eq.  (A-2), 


hu,3  =  0,  9  molecule. 


(A-5) 


To  obtain  the  four  enthalpies  for  the  mixture,  a  mass  average  is  performed: 


h,  =  ^— q  =  i,r.v,€. 
,=i  P 


Then,  the  mixture  enthalpy  is  given  by; 


A=y;^'... 
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The  energies  and  their  associated  temperatures  are  related  by  specific  heats  according  to: 


q  =  r.i'.p,  s  =  1 . ns. 


Also,  the  enthalpies  and  the  temperatures  are  related  by; 


q  =  t,r,v,6,  s  =  1 . us. 


Here,  ^  and  ^  are  specific  heats  per  unit  mass  at  constant  volume  and  pressure,  respectively. 
From  statistical  mechanical  considerations  it  can  be  shown  that  (Ref.  32) 


r'*  —  ^  ^  t2  ^ 


MsdT,,,  T"-*  dT,,s  J’ 


q  =  t,r,  p,€,  s  =  1, . . . 


(A-6) 


Then,  by  Eq.  <A-3),  it  follows  that: 


f  cr,  <i  i 
s  -  J  n 

P.tf  1  r"  +  —  q  =  t 

[  " «  .VI, 


(A-7) 


As  with  the  enei^es  and  enthalpies,  the  specific  heats,  (fy  and  C^,  for  a  given  species,  s,  are 
assumed  to  be  the  sum  of  translational,  rotational,  vibrational,  and  electronic  parts: 


ft  _  1  fa  A.  pa  A.  pa 

'^  L'  —  e.i  I  ^  i',r  T  -1- 

pa  _  p»  I  r-s  I  pa  I  pa 

p.<  p.r  '  p,v  '  p,e> 


S  =  1,  .  . .  ,US. 


(A-8) 


To  obtain  the  mixture  specific  heats  for  the  different  modes,  a  mass  average  is  performed: 


p  _  ?±p^ 

8=1  “ 


Pa  ^8 
8=1  r 


q  = 


Then,  the  mixture  specific  heats  are  given  by: 


/IB  _ 

c,  =  Y.^ci 
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In  complete  thermal  equilibrium,  when  ail  temperatures  are  equal,  say  T,.,  =  T,  Eqs.  (A-6)  and  (A-7) 
can  be  evaluated  explicitly  to  give  the  species  specific  heals  according  to  Eq.  (A-8).  The  following 
interpolation  is  used  to  approximate  this  calculation  as  a  function  of  the  equilibrium  temperature; 


'  i=l 


C’lT)  -  —Ya  T-^-~ 


» i=i 


(A-9) 


Also,  in  certain  cases  where  the  temperatures  are  distinct,  explicit  calculations  can  be  performed 
easily  to  obtain  the  specific  heats.  For  example,  the  translational  and  rotational  specific  heats  are 
independent  of  temperature: 


311 

2M, ' 

■  2M,' 

(A-10) 

ra  _ 

■  2M, 

=  C*  . 

p,r' 

(A-ll) 

Here,  r^is  the  number  of  rotational  degrees  of  freedom  for  a  species,  s.  In  particular,  for  diatomic 
molecules,  2.  Next,  though  it  requires  more  effort,  the  electronic  specific  heats  can  be  calculated 
from  Eq.  (A-6)  after  tabulating  a  sufficient  number,  say  /,,  of  energy  levels,  ,  and  degeneracies. 

The  result  is; 


(A-12) 


Finally,  because  of  the  complexity  of  calculating  the  vibrational  specific  heats  by  Eq.  (A-6 )  for  real 
molecules,  the  following  approximation  is  used: 


/■•a 


=  -  CtJT,.,)  -  -  C’lAT,.,) 


c 


<0 

p.v 


CUT,,a) 


(3+  rjT? 
2Ms 


c 


•s 


(A-13) 


Here,  the  nonvibrational  parts  are  subtracted  away  from  the  total  specific  heat  for  the  species,  after 
evaluating  all  these  specific  heats  at  the  vibrational  temperature  in  Eqs.  (  A-9)  -  (A-12). 


Since  the  free  electron  has  only  translational  energy  or  enthalpy,  its  specific  heats  have  only 
a  translational  part.  Thus,  with  s  =  e,  Eq.  (A-8)  becomes: 
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/'<f~  _  _ 

— 


dU 

2M/ 


n^~  _  /'>e  _ 

p  -  '~p,f  - 


hll 

2M/ 


(A-14) 


For  the  temperatures,  q  =  t,r,  v,  e,  there  are  associated  thermal  conductivities,  q 
=  t,  r,  V,  e,  respectively.  First,  the  translational  thermal  conductivities  are  given  by 


«<.*  = 


15 


7^ 


(A-15) 


Here,  a,  is  the  collision  diameter  for  the  pure  species,  s,  and  Q  is  the  collision  integral  defined  as 
the  weighted  average  of  a  collision  cross  section  of  the  form 

-  cos' 

r’  »  - ; - ’ 

f  I  e”"^  7^'  ■’■^11  —  cos'  v)sin  \d\(h 

Jo  Jo 


_  [  MM 


<A-16) 


Here,  is  the  collision  cross  section  for  the  s-r  collision  pair,  %  is  the  scattering  angle  in  the  center 
of  mass  system,  g  is  the  relative  velocity  of  the  colliding  particles,  and  y  is  the  reduced  velocity.  The 
expression  in  Eq.  (A-1 5)  is  approximated  within  NEQPAK  (Ref.  23)  by  the  following  interpolation: 


f*1.3  — 


1.5  71 


4 

The  other  thermal  conductivities  are  given  by: 


exp 


s 

r 

Li=l 


^  6,-,,  (In  r, 


(A-17) 


''9, a  — 


q  =  r,v,e. 


(A-18) 


Here, 


Vs  =  Vss{Tt,s) 


is  the  self-diffusion  coefficient,  where,  in  general,  the  binary  diffusion  coefficient  is  given  in  terras 
of  a  temperature,  T,  as: 


'P.,(ri 


3^-r^  iTliMs  +  ■Vf,.T 

SpSl[YhT)\  'ivMsMr 


(A-19) 
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This  expression  is  approximated  within  NEQPAK  (Ref.  23)  by  the  interpolation. 


S  (T) 

V„(T)  =  10.1325-2^  exp 


j-i 


4=1 


where  S^(T)  is  the  electrostatic  shielding  ratio  (Ref.  27), 

1  s  and  r  neutral 


S,AT]  =  { 


1 


In  A 


K*) 


«  or  r  charged 


(A-20) 


and  the  shielding  function  A(T)  is  defined  as 


A(r)  = 


2.09  /'T'* 


fT‘*\  1.52  /T^ 

VPe-/  ^  10^®  \Pc- 


and  pg.  is  the  partial  pressure  of  the  electron  in  atm.  In  case  a  collection  of  species  has  the  same 
temperature  for  a  given  energy  mode,  an  associated  mixture  thermal  conductivity  can  be  defined  in 
terms  of  thermal  conductivities  for  the  pure  species.  Specifically,  suppose  all  species  but  the  free 
electron  have  the  same  translational  temperature,  T.  Hien,  the  translational  thermal  conductivity  for 
this  mixture  is  given  by: 

■ill  Ll,ris  \i 


— 


det 


det 


•in«.l 

A'l 


0 


3,r?£*  ~ 


Here  the  terms  depend  on  the  pure  species  thermal  conductivities,  but  the  explicit  form  is  given 
below  only  for  the  diagonal  terms  to  fully  specify  the  approximation. 


where  Xj  is  the  mole  fraction  of  species,  s.  The  diagonal  terms  are  given  by 

VAr  {^M]  +  ^M't  +  4MrM,A,r) 


(A-21) 


y/2{Mr  +  A, 


Mj  -VI J 

Z-F  A-  B 

^  ■*  *r  T  ^  Da 


1  i 


'i.r 
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where  in  generaJ, 


{T) 


F„  = 


\ 


0.2.2) 


(r) 


B,r 


ni^ 


m 

{T) 


Following  the  method  of  Armaly  and  Sutton  (Ref.  29),  these  terms  are  approximated  as  follows: 

'  0.18  s  =  H,r  = 

_  I  _  J  0025  ^  =  H(,r  = 

A»,-  -  Ar,  -  <  ^  ^  ^  r  =  »-^  ' 

1.2-5  otherwise  (A-22) 

and 

F,r  =  Frs  =  h  (A-23) 


0.2  41  =  atom  or  molecule,  r  =  e~ 

0.15  s  =  atom  or  molecule, r  =  ioa 
0.78  41,  r  =  atom  or  molecule 

1.0  otherwise  (A-24) 


Thus,  Eq.  (A-21)  can  be  viewed  as  giving  a  mixture  thermal  conductivity  according  to  a  mixture 
function,. Vas>  due  to  Armaly  and  Sutton,  that  depends  on  pure  species  thermal  conductivities. 
In  particular,  if  the  translational  temperatures  for  all  species,  except  the  free  electron,  are  the 
same,  then 

fit  —  A  As(  •  •  -  •  >■■•)•  (A-25) 

In  a  similar  manner,  other  mixture  thermal  conductivities  can  be  defined,  provided  the  constituents 
have  the  same  temperature  for  the  corresponding  energy  mode.  Specifically,  if  the  rotational 


temperatures  are  the  same  (there  is  none  for  the  electron), 

fir  =  Aa,?( - - )‘  (A-26) 

If  the  vibrational  temperatures  of  the  molecules  are  the  same, 

K,j  —  Aas(  .  • .  •  4>i:,9=iiu)lM»Ie*  •  •  •  )•  (A-27) 
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If  the  electronic  temperatures  are  the  same, 

- ).  (A-28) 


Since  the  free  electron  has  only  a  translational  temperature,  i.e.,  7^  =  Tf  its  thermal  conductivity 
is  given  by  Eq.  (A-IS), 


1.5  n 
4 


exp 


^6,,-(lnr,_)*-* 


t=l 


(A-29) 


Finally,  if  all  temperatures  are  the  same  and  equal  to  T,  the  species  thermal  conductivities  are 
determined  by 


1.5  n 


S— #■ 


1=] 


(A-30) 


Here,  the  first  expression  shows  the  translational  part.  TTie  second  expression  shows  the 
translational  part  of  the  specific  heat  subtracted  away  so  that  the  expression  gives  the 
nontranslational  part  of  the  thermal  conductivity.  Then,  the  mixture  thermal  conductivity  is  given, 
according  to  the  method  of  Aimaly  and  Sutton,  as: 


n 


(A-3I) 


A  similar  procedure,  due  to  Armaly  and  Sutton  (Ref.  28),  is  used  to  determine  the  mixture 


viscosity,  p.  First,  the  pure  species  viscosity,  for  a  species,  s,  is  defined  by 


15  Tv 


(A-32) 


where  K,,  is  defined  in  Eq.  (A-15).  Thus,  by  Eq.  (  A-17),  p,  is  approximated  by  the  interpolation, 

5 


1=1 


(A-33} 


Then,  the  viscosity  for  the  mixture  is  given  by 


■  -ffll  • 

1 

Hns,na 

\,n« 

>^1  ■ 

■  \na 

0 
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Here  the  terms  H„  depend  on  the  pure  species  viscosities,  but  the  explicit  form  is  given  below  only 
for  the  diagonal  terms  to  fully  specify  the  approximation, 


The  diagonal  terms  are  given  by 


n*  .  3 

.=  1 


(A-34) 


//«  =  ^  +  E 

it 


Vs  Vr 


^1.  2>/2(,Vf.  +  -V(,)§ 

where  and  are  given  by  Eqs,  (A-22)  -  (A- 24). 


Mr] 

1 

t 

1 

Ms\ 

+ 

_ 1 

D 

■i2 


In  the  remainder  of  this  report,  it  is  assumed  that  the  translational  and  rotational  energies, 
and  e„,  s  *  e‘,  are  distributed  with  the  same  temperature,  T,  i.e., 

a  ^  f  . 


Also,  the  electronic  excitation  energies,  s*  e,  are  assumed  to  be  distributed  with  the  same 
temperature,  T„  i.e., 

Tf,f-Ti.  s^e~. 

Recalling  7”,.=  F,,.,  and  using  Eq.  (A-3 ), 


K- 


=  ft- 


+ 


nr,- 
M,-  • 


Since  there  is  a  single  translational  temperature,  Eq.  (A-3)  shows  that  the  translational  enthalpy  and 
energy  for  other  species  are  related  according  to: 

UT 

hus  =  et.,  +  ^.  (A-35) 

Al, 

With  the  above  definitions,  the  conservation  equations  can  now  be  given.  First,  the 
conservation  of  mass  for  each  species  is  expressed  as: 


0 


j  =  1, ....  nil 


(A-36) 
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where  is  the  diffusion  velocity  of  a  species,  s.  Also,  a>,  is  a  source  term  that  depends  on  the 
partial  densities  and  the  translational  temperature  in  such  a  way  that  it  vanishes  when  equilibrium 
conditions  are  achieved.  Specifically,  the  source  terms  are  constructed  as  follows.  The  reactions  on 
which  they  are  based  can  be  written  in  the  form 

nj  ns 

n=l  jn=l  (A-37) 


Here,  and  are  stoichiometric  coefficients  for  the  rth  reaction,  M, 

represents  one  mole  of  species  s,  and  nr  is  the  number  of  reactions.  The  rate  of  disappearance  of  a 
species  s  due  to  the  rth  reaction  is 


■1  =  1 


11=1 


(A-38) 


while  the  production  rate  of  a  species  s  due  to  the  rth  reaction  is 

n»  114 

Sr.i  =  JJ  Om  *  +  l^r.a^'r  JJ  7 


m=l 


The  net  rate  of  change  of  species  s  due  to  all  reactions  is 

7ir 

<^s  =  M» 

r=l 


(A-39) 


(A-40) 


Like  most  aerothermochemical  models,  NEQPAK  assumes  the  modified  Arrhenius  form  (Ref.  32) 
for  chemical  reaction  rates 


=  ArT^^e:rj}{-C,./T),  (A-41) 

where  T  is  the  translational  temperature  of  the  gas.  The  rates  defined  in  Eq.  (A-41)  were  obtained 
from  fits  to  experimental  data  that  were  collected  under  conditions  of  thermal  equilibrium. 
However,  when  the  reaction  rate  under  question  represents,  for  instance,  a  two-body  dissociation 
reaction,  then  the  reaction  rate  also  depends  on  the  vibrational  temperature  of  the  gas.  One  of  the 
mote  popular  approaches  used  to  account  for  thermal  nonequilibrium  is  that  of  Park  (Ref.  33),  who 
replaced  the  translational  temperature  in  Eq.  (A-41)  with  a  generic  temperature  or  average 
temperature,  T^.  Park  defined  7,  by  the  relation 

Tf,  =  TyT^  ^;0  <  O’  <  1,  (A-42) 

where  a  was  chosen  to  reproduce  experimental  data. 
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The  reverse  reaction  rate  coefficient  is  related  to  the  forward  rate  by 


SFr 


where 

A’;’  = 

is  the  equilibrium  constant  for  the  rth  reaction  and 


B=i 

718 

SFr  = 


9^\ 


(A.43) 


(A-44) 


(A-45) 


In  Eq.  (A-44),  =  1 .01325  x  10^  Pa.  Also  C®  is  the  Gibbs  free  energy  for  species  s  at  the  given 

temperature  and  a  pressure  of  I  atm.  The  Gibbs  free  energies  are  computed  from  curve  fits  taken 
from  Ref.  27.  This  interpolation  takes  the  form, 


6'J{T)  =  7?.r 


r*  -I-  -  ttr,. 


Next,  the  conservation  of  mixture  momentum  is  expressed  as: 


pUfllj  +  iij 


t  =  1,2,3 


(A-46) 


Here,  the  mixture  viscosity,  p,  is  given  in  Eq,  {  A-34  ). 

The  vibrational  energy  associated  with  each  molecule  is  conserved  according  to: 


0  ,  ^ 
01 


Oxj 


fl  =  molecule. 


(A-47) 
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where 


—  pa' 


T’t-u.j 


+  iOgDa 


Here,  e*  ^  and  are  the  vibrationa]  energies  of  species  s  evaluated  at  the  translational 
temperature,  7,  and  the  electronic  excitation  temperature,  T„  respectively.  Also,  and  are 
the  translational-vibrational  energy  relaxation  time  and  electronic-vibrational  energy  relaxation 
time,  respectively,  for  molecular  species  s.  Finally,  t>e  is  the  average  vibrational  energy  per  unit 
mass  of  molecule  s,  which  is  created  at  rate  09,. 


Since  all  the  electronic  excitation  energies,  are  assumed  to  be  distributed  with  the  same 
temperature,  7,,  the  conservation  of  the  mixture  electronic  excitation  energy,  e„  can  be  expressed 
in  terms  of  this  single  temperature  as  follows: 


£ 

dt 


0 


i=i 


—  Qratl 


(A-48) 


where  is  given  in  Eq.  (A-28).  Also,  is  the  radiative  energy  transfer  rate.  Since  the 
translational  temperature,  7^,.  =  7,.,  for  the  free  electron  is  assumed  to  be  different  from  that  for  oth¬ 
er  species,  the  conservation  of  electron  energy  must  be  written  separately  as; 


U  L 


p^-e^~Uj 


where 


j=i  ■> 


^,-  =  3p,-n{T-T,-]  Y, 


a^e- 


M 


-  -  H  -  Y.  P‘ 

’  »=ioil  S3:ilIol. 


r._ 


(A-49> 


and  1^.  is  given  in  Eq.  (A-29).  Also,  is  the  effective  collision  frequency  for  electrons  and  heavy 
particles  in  electronic-translational  energy  relaxation.  Then,  is  the  molar  rate  of  production 
of  species  s  per  unit  volume  by  electron  impact  ionization,  and  /,  is  the  fust  ionization  energy  of 
species  s  per  kg-mole. 


Finally,  define  the  total  energy  as  the  sum  of  internal  and  kinetic  energy; 

3 

=  e  +  0.5  Ilf. 

IBl 
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This  is  conserved  according  to  the  following; 


0  4^  d 

^(/>eT}  +  E7»7' 


L 


dT 


(/)CT  +  phi}  -  (fit  +  -  2^  Kf.5 

^  i=inol. 


fir...  .  aTe  .  OT,- 

?v  ■  .V  ' 


Ji=l 


d-Vj  'OXj  *■  C>.Tj 

/  duj  duj  ^  2  ^ 


\  ^  tfsion  / 

where  tc^  and  aie  given  in  Eqs.  (A-25)  and  (  A-26). 


(A-50) 


Now,  assume  that  the  electron  energy,  and  the  mixture  electronic  excitation  energy,  e,, 
are  distributed  with  the  same  temperature,  T^,  i.e., 


Te  =  T,-  =  r. 


Therefore,  the  conservation  of  the  two  energies  can  be  expressed  in  a  single  equation  in  terms  of  a 
single  temperature.  For  this,  define  the  electron  and  electronic  excitation  energy, 

f£  =  e, 

P 

Also,  certain  molecules  may  be  assumed  to  have  the  same  vibrational  temperature.  These  will  be 
called  type  I  molecules,  and  their  vibrational  temperatures  are  assumed  to  be  equal  to  the 
translational  temperature,  i.e., 

T  =  r,.,*,  s  =  type  1  molecule. 


The  remaining  molecules  with  distinct  vibrational  temperatures  will  be  called  type  II  molecules. 
Finally,  the  degree  of  charge  separation, 

j=ion 

is  assumed  negligible  and  is  set  to  zero  in  the  appropriate  conservation  equations. 


Now  the  conservation  equations  are  given  as  follows.  First,  the  conservation  of  mass  for  each 
species  is  expressed  as 


0 


d 


j=i  •' 


n 

P.“J  - 


=  CJ, 


9~  1. . . . 


(A-51) 


89 


AEDC-TR-94-18 


Note  that  in  this  equation  and  in  ones  that  follow,  the  diffusion  velocity  of  species  s  is 
approximated  as 


^  P3  Olj 


(A-52> 


where  X  is  the  mole  fraction  of  species,  s,  and  £),  is  an  effective  diffusion  coefficient  given  by 

In  genera],  the  diffusion  velocity  for  a  given  species  depends  on  concentration  gradients  of  other 
species  as: 


VluT 


J 

where  Dg  is  the  thermal  diffusion  coefficient  for  species  s.  This  dependence  does  not  appear  in  the 
approximation  of  Eq.  (A-S2  ). 


Next,  since  the  number  density  of  the  free  electron  is  assumed  approximately  equal  to  the 
total  number  density  of  ions,  the  conservation  of  mixture  momentum  is  expressed  as: 

ii  0  i.  (  '^  ^  / On,  df/Al 

^ 5  557  + 'H"  +  3"  £  ^  j  ■ "  ^  ^  j  J " " 

/  =  1.2.3  (A-54) 

The  vibrational  energy  associated  with  each  type  II  molecule  is  conserved  according  to: 


d  ^  ^  d 

Of  L' 


,  n 


=  U) 


where 


dxj  '""'’■'‘'‘cl.r, 

3  =  type  II  molecule, 


(A-55) 


~  Pa' 


^v.a  ~ 


Pa- 


^v,a  ~  ^v,a 


Te-i 


Since  the  mixture  electronic  excitation  energy,  e^,  and  the  electron  energy,  e,.,  are  assumed 
to  be  distributed  with  the  same  temperature,  Te,  the  conservation  of  the  electron  and  electronic 
excitation  energy,  e^,  can  be  expressed  in  terms  of  this  single  temperature  as  follows; 
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;=1  ^ 


-  (*>'e  +  Kf-  )~Q^  -  py'.hs.sDa^^ 


3=1 


d.r 


where 


and 


Oil: 

=  -l^Pr-^  +  u;£  -  Qr«| 

j=l 


(A-56) 


f^E.. 


=  / 

^  ^le.s  C 


_  s  =  e 
otherwise 


^E  =  ip,-n(T-TE)  E  ^  -  E  W.-  E  ^ 

,/p-  «i«i  - . ' 


#•*  -r 

^  i',<  '  »'- 


jsmul. 


Finally,  the  total  enei^y  is  conserved  according  to  the  following: 
i)  0  ,  ,  f  ,  ,  v~‘  \  V'  OTi,_i 

57i'’'Ti  +  Zar  h  +  “’+  L  L  ''■■■■ife- 

/=!  J  L  '  ^snw*  I  '  4=iiiol.  II  ^ 

OTe  y  _  0\$  ( Oui  2  ,  r— ^  « 

-(K.  +  Ke-  )7^  -  P'E  S  ^  ~ 

(A-57) 

Now,  assume  that  the  vibrational  energies  of  all  molecules,  s  =  molecule,  are  distributed 

with  the  same  temperature,  T„  i.e., 

Tu3  =  Tu,  s  =  molecule. 


Therefore,  the  conservation  of  vibrational  energies  can  be  expressed  in  a  single  equation  involving 
the  mixture  vibrational  energy,  e„  and  the  single  temperature,  Ty. 


Now  the  conservation  equations  are  given  as  follows.  First,  the  conservation  of  mass  and  of 
mixture  momentum  are  shown  in  Eqs.  (A-54)  and  (A-36).  The  mixture  vibrational  energy  is 
conserved  according  to: 


0  ^  0 
j=l  *' 


OT,  ^  ,  _  0\s 


(A-58) 


where  K,,  is  given  in  Eq.  (A-27),  and 


lu* 


E 

«=11U>1. 
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The  conservation  of  die  electron  and  electronic  excitation  energy  is  expressed  as: 


(9,  ,  ^  0 

j=i  ■> 


£vp  IIJI  ft 

pesUj  -  (Kg  + 
dUj 


(A-59) 


where 


LUUJ  ^ 

Pt~  .  "!■<*''£  ~  Qrad 
j=l  J 


wj:  =  V,-Kir-Ii:)i:  ft..,/.-  5; 

»^€-  ^  «=jon  j=njcil.  c-(SA 

Finally*  the  total  energy  is  conserved  according  to  the  following: 


v  \  % 


0 


dT  dT  fyn  ■) 

(/>fT  -  (Kj  +  - -  {Hf  +  Ke-)-T~  -  /’Y' 

dxj  Oxj  Oxj  '  ^  Oxj 


f  dui  Ouj\  2  ^  ^ 

j  +  s'"*'*'-''  g  airj  = 


(A-60) 


Now,  assume  that  the  electron  and  electronic  excitation  energy,  and  the  mixture 
vibrational  energy,  e„  are  distributed  with  the  same  temperature,  Ty,  i.e., 

Te  =  T,  =  Ti'. 


Therefore,  the  conservation  of  the  two  energies  can  be  expressed  in  a  single  equation  involving  the 
mixture  vibrational  and  electronic  energy, 

ev  =  eE  +  fi-. 


and  the  single  temperature,  Ty. 


Now  the  conservation  equations  are  given  as  follows.  First,  the  conservation  of  mass  and  of 
mixture  momentum  are  shown  in  Eqs.  (A-S4)  and  (A-36).  The  mixture  vibrational  and  electronic 
energy  is  conserved  according  to: 


J=1  J 


.dTv 


ftevllj  -  (Kf  +  Kf  +  K,-  -p'^ 


dxj 


a=l 


.;=! 


£i!i 

Oxj 


"I"  “  Qrad 


(A.61) 
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where 


U'.* 


{  lh,s  +  h,.,» 

[ 


s  =  c~ 

»  =  inol^rule 
otherwise 


and 


\jjy  —  tjJv  "h  '^E 


Finally,  the  total  energy  is  conserved  according  to  the  following: 


Now,  assume  that  all  energies  are  distributed  with  the  same  temperature.  The  conservation 
equations  are  given  as  follows.  First,  the  conservation  of  mass  for  each  species  is  expressed  in  Eq. 
(A-36).  The  conservation  of  mixture  momentum  is  expressed  in  Eq.  (A-54).  Finally,  the  total 
energy  is  conserved  according  to  the  following: 


(A-63) 


Here,  tc  is  given  in  Eq.  (A-3 1 ). 
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APPENDIX  B 
CELL  GEOMETRY 

The  method  by  which  the  flow  field  is  discretized  into  a  three-dimensional  computational 
domain  is  described  in  this  appendix.  The  numerical  algorithm  employed  by  NEDANA  is  a  finite- 
volume  scheme.  Traditional  finite-volume  schemes  are  constructed  on  computational  grids  where 
the  dependent  variables  (Q)  are  stored  at  cell  centers,  and  the  independent  variables  (jc,  y,  z)  are 
known  at  cell  comers.  Employing  the  notation  of  Gnoffo  (Ref.  21),  (J,  K,  L)  denote  the  indices  of 
the  cell  centers,  and  (/,  k,  1)  denote  the  indices  of  the  cell  faces.  The  two  indexing  systems  are  related 
by  (/,  K,  L}  =  (j  +  1/2,  k  +  1/2,  /  +  1/2).  In  traditional  schemes,  the  cell  comers,  (j,  k,  [),  form  the 
vertices  of  a  computational  cell,  and  the  dependent  variables  represent  the  average  of  the  variables 
over  the  cell  volume  (7,  K,  L).  One  disadvantage  of  this  grid  technique  is  that  the  grid  must  be 
shifted  to  output  the  (jc,  y,  z)  location  of  the  dependent  variables.  Because  NEDANA  interacts  with 
software  that  requires  the  dependent  and  independent  variables  to  be  written  at  the  same  location, 
an  alternative  grid  system  was  adopted.  In  this  system,  the  computation  grid  is  constmcted  from  the 
cell-centered  coordinate  system.  In  effect,  the  coordinate  system  is  shifted  by  half  an  index.  The 
areas  and  volumes  are  calculated  on  the  (7  -  1/2,  IIT  - 1/2,  L  - 1  /2)  coordinate  system  and  stored  before 
the  flow  solver  is  executed.  This  shift  is  made  possible  because  the  location  of  the  cell  comers  is 
never  explicitly  used  in  the  finite- volume  scheme.  As  long  as  the  cell  volumes  and  areas  are 
calculated  correctly  with  respect  to  the  celt  centers,  the  coordinate  system  chosen  is  irrelevant.  This 
shifting  is  transparent  to  the  user  of  NEDANA  because  all  of  the  cell  geometry  logic  is  performed 
by  a  preprocessor  code. 

The  shifting  does,  however,  affect  the  treatment  of  the  boundary  conditions.  Unlike 
traditional  finite-volume  schemes  where  domain  boundaries  correspond  to  cell  faces,  the 
NEDANA  boundaries  correspond  to  cell  centers.  Therefore,  the  boundary  conditions  can  be  set 
directly.  Figure  B-1  is  a  schematic  of  the  NEDANA  grid  system  for  a  single  plane.  Note  that  at  the 
boundaries,  the  cell  face  areas  become  cell-centered  areas.  This  psuedo-area  is  appropriate 
because  no  flux  evaluations  are  performed  at  the  boundaries.  The  areas  are  only  used  for 
calculating  normals  to  the  surface  and  setting  boundary  conditions.  In  Fig.  B-2,  the  geometry  of  a 
single  computational  cell  is  shown.  In  this  figure,  the  two  coordinate  systems  and  the  relative 
positions  of  the  cell  faces  and  the  cell  center  can  be  more  easily  seen.  Note  that  the  cell  faces  lag 
the  cell  centers  by  half  an  index. 
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The  calculation  of  the  cell  face  areas  and  volumes  is  performed  entirely  by  the  NEDANA 
preprocessor  code.  The  methodology  employed  follows  directly  from  that  of  Gnoffo  (Ref.  21). 
First,  define  the  location  of  the  cell  comers  with  position  vectors  in  the  Cartesian  coordinate 
system.  Let 


where 


dfi 

= 

df2 

= 

= 

df4 

= 

dr’s 

= 

dfe 

= 

dfj 

= 

(B-l) 

= 

+  yiy  +  • 

(B-2) 

The  directed  cell  areas,  a,  and  the  cell  volume,  V,  can  now  be  defined  by 


^  dfi  X  dr2 

-  2 

_  dfa  X  dfi 

=  2 

,  dr's  X  dfa 

-  2 

(B-3) 

(B-4) 
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The  projected  area,  o,  with  magnitude  a,  is  normal  to  the  cell  face  and  points  in  the  direction  of 
increasing  J,  K,  L,  respectively.  In  the  NEDANA  notation,  the  cell  face  areas  are  written 

+  <r^iz 

^vi.K.L  =  er^ix  +  +  cr=l^ 

^Cj.A'.L  = 

where  the  terms  (f,  c^,  and  <f  denote  the  magnitude  of  the  area  as  projected  onto  the  Cartesian 
coordinates.  The  transformation  metrics  in  a  finite-volume  scheme,  such  as  ^  ^2,  are  expressed 

as  the  ratio  of  cell  face  areas  to  cell  volumes.  Equation  (B-4)  is  first-order  accurate  with  respect  to 
a  cell  face  or  a  cell  volume.  Gnoffo  states  that  this  formulation  introduces  errors  when  crossing  an 
axis  singularity.  Gnoffo  proposes  a  second-order  accurate  expression  that  uses  symmetric  averages 
of  differences  about  the  cell  center.  The  volume  can  now  be  expressed  as 


+  ^04-1, A-.L  + 

^  (B-6) 


where 


lifs  =  [(»^  +  x,,  +  4-  +  ^ij  + 


(B-7) 


and 


i^dj,A',L  - 

2 

fi  1  — 

[*'j]j+l,A',L  _ 

rs-1  — 

2 

(B-8) 


where  the  generic  differences  r  are  second-order  accurate  with  respect  to  the  cell  centers,  and  the 
dummy  variable  s  represents  the  independent  variables  x,  y,  z.  The  s  differences  are  second-order 
accurate  with  respect  to  the  cell  faces  and  have  the  form 


[''‘dj.A'.L 


i 

2l^J+l.t,l  -  Sj.kJ  +  Sj+i,k+l.l  -  9j^li+lj) 
2(*J.t'+l.l  ~ 

+  "V*+1J+1  “  ^j.fr+1./) 
^(fj.U+l  -  »j,U  +  »j+\,kJ+Ti  -  8j+l.t./}. 


(B-9) 
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The  transformation  metrics  are  defined  by 


& 

Vx 

Cx 

Vi 

Vv 

— 

Vv 

V: 

. 

vc 

’<  . 

(B-10) 


In  finite-volume  notation,  these  terms  are  written  as  the  ratio  of  cell  face  areas  to  cell  volumes 


V4j,A',L  = 

^0,K,l  = 


^J,K.L 

^vJ.k,L 

Vj.A.Z, 


(B-11) 


Equations  (6-1 1)  are  only  first-order  accurate.  The  reason  being  that  a  is  second-order  accurate 
with  respect  to  cell  faces  but  only  first-order  accurate  with  respect  to  cell  centers;  and  conversely. 
V  is  second-order  accurate  with  respect  to  the  cell  centers  but  only  first-order  accurate  with  respect 
to  cell  faces.  Gnoffo  suggests  using  a  symmetric  average  of  cell  faces  and  volumes  to  obtain 
second-order  accurate  metrics  for  the  evaluation  of  viscous  terms.  These  second-order  accurate 
metrics  take  the  form 


^V},K.L 


vo./f.t 


4  Vj.a'.jlVj-i.a.l 

+  ^tj+l.A-l,L)  +  V  J.A'-l, //(<?!£  j,  A, L  + 

Vj,A',l{<T(j,A',L-1  +  +  ^ij+l,h\L) 

4  Vj,a.lVj,a.l-i 
4  Vj,A',lVj_i,a.L 

Vj,A,A(gt)J.l:-l.L  +  ^JiJ,k,L)  +  ^J,h'-l,L{^r,J,k,L  + 

4  Vj,a',lVaA-1.L 

Vj,A.A(gt)  AA-.L-l  +  +  ^r,J,k+l,L) 

4 

J-l.AM  +  -I-  V./-1,A',a(^CJ,A,/  + 

4  Vj,A'.LVy_i,A',L 

J-l.A./  + 

4  Vj.a-,lVj,k-i.l 

^J,K,l{^CJ.KJ~^  +  +  Vj.A',L-l(^CJ,AV  + 

4  Vj.axVaa,l-i  (B-12) 
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APPENDIX  C 

NUMERICAL  BACKGROUND 

The  purpose  of  this  appendix  is  to  discuss  basic  issues  crucial  to  the  selection  of  a  finite- 
volume  scheme  for  the  solution  of  the  nonequilibrium  flow  equations  shown  in  Eq.  (43).  The  topics 
covered  include  first,  the  properties  of  schemes  that  guarantee  convergence  of  approximate 
solutions:  conservativeness,  consistency,  and  total  variation  stability.  Then  the  flux-limited  and  the 
slope-limited  approaches  to  constructing  high-resolution  convergent  schemes  are  introduced.  The 
examples  given  to  illustrate  these  two  approaches  are  Yee’ s  symmetric  TVD  scheme  and  van  Leer’s 
MUSCL  scheme,  respectively.  However,  a  primary  objective  of  the  work  documented  in  this  report 
is  that  the  numerical  method  selected  be  not  only  fast,  accurate,  and  stable,  but  also  easy  to  under¬ 
stand  and  modify.  With  this  in  mind,  Jameson’s  flux-limited  dissipation  model  is  developed  as  an 
alternative  to  the  previously  defined  techniques.  Most  of  this  material  is  presented  for  scalar  model 
problems  but  also  is  briefly  generalized  to  a  vector  setting.  For  a  more  detailed  discussion  of  the 
issues  addressed  in  this  appendix,  see  LeVeque  (Ref.  61). 

For  simplicity,  first  consider  methods  for  solving  the  model  problem, 

Qi  +  fx  =  0.  (C-1) 

The  task  is  to  approximate  the  solution  to  this  problem  on  a  space-time  grid,  {(x^,  i")},  with  grid 
values  [^ "}  such  that  q"  «  q{xj,  ^).  Although  all  grid  values  are  known  only  at  grid  points,  the 
definition  of  certain  numerical  constructs  will  be  motivated  by  focusing  on  activity  at  the  midpoints, 
Xj+m  =  l/2(x^  +  ^^i)-  These  will  also  be  referred  to  as  the  interface  points  between  adjacent  spatial 
cells.  The  jth  such  cell  is  defined  as  the  set  of  points  for  which  x  is  between  Xj.,a  and  Xj^,^. 

It  will  also  be  necessary  to  consider  the  performance  of  approximation  schemes  on  space-time 
grids  of  ever  increasing  refinement.  Specifically,  it  is  crucial  that  the  grid  values  {qj^}  approximate  a 
solution,  q,  with  arbitrarily  small  error,  provided  Ar  =  and  Ar  =  Xj+iq-  Xj.ia  are  sufficiently  small. 

This  convergence  is  theoretically  guaranteed  of  finite-difference  schemes  that: 

1 .  can  be  written  in  conservation  form, 

2.  are  consistent,  and 

3.  are  total  variation  diminishing  (TVD). 

These  three  conditions  are  explained  in  detail  below. 
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C-1.0  CONSERVATION  FORM  AND  NUMERICAL  FLUX  FUNCTIONS 


Before  deriving  a  finite-difference  approximation  to  Eq.  (C-l),  note  that  a  physically  rele¬ 
vant  solution,  q,  to  Eq.  (C-I)  may  develop  discontinuities  even  in  cases  in  which  the  initial  data, 
q{x,  0),  are  smooth  (Ref.  62).  Such  a  q  would  not  satisfy  the  differential  equation  in  a  strict  sense. 
Thus,  it  is  productive  to  express  Eq.  (C-l)  in  a  weaker  form,  where  no  differentiability  is  required 
of  the  solution.  This  is  accomplished  by  integrating  the  differential  equation  over  an  arbitrary  space- 
time  cell,  say  [jCj,afj]  x  {/j,  /j].  Specifically,  integrating  explicitly  with  respect  to  t,  and^ explicitly 
with  respect  to  x  gives  the  integral  equation; 

/  -  q{x,U)]dx  +  f  [/(9(a:2,t))  - =  0- 

Jxi  ■'tl 

Notice  that  no  derivatives  appear  in  this  equation.  In  particular,  the  meaning  is  clear  even  if  a 
solution  is  discontinuous.  A  function,  q,  that  satisfies  this  equation  for  all  jr„  Xj,  ti  and  is  said  to 
be  a  weak  solution  to  Eq.  (C- 1 ).  Toward  a  finite-difference  approximation  to  this,  let  the  integration 
cell  be  [jCj.1/1.  x  [r,  so  that  it  is  aligned  with  a  space-time  grid  cell.  Then  the  above 
equation  can  be  written  in  the  finite-volume  form: 

[,-«+!  _  ^^]Ax  -h  =  0. 

Here,  ^  and  q"  * '  denote  averages  of  the  solution  over  the  spatial  interval,  at  times 

t  =  f  and  /  =  r*‘,  respectively.  Also,^.,n  andjJ^Q  denote  averages  of  the  flux  over  the  temporal 
interval,  (f,  at  cell  interfaces  x  ~  Xy.1/2  and  x  =  respectively.  This  calculation  motivates 
the  distinction  of  finite-difference  schemes  in  the  following  conservation  form: 


-17 1 

At  Ax 


=  0. 


(C-2) 


Compare  the  last  two  equatmns.  Here,  the  grid  value,  q”,  approximates  the  spatial  average,  q”. 
Also,  the  temporal  average,  />4.iy2  is  approximated  by  the  numericaf^ux/Hncrion,  Asample 
/  is  given  by: 


(C-3) 


which  leads  to  the  central  difference  approximation  to/^’ 


Ax 


‘2Ax 
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In  general, 1/2  ^  function  of  a  certain  number  of  grid  values  on  opposite  sides  of  the 

interface  at  Also,  these  could  be  taken  from  time  level  f*  or  r‘*‘.  In  the  example  shown  in 

Eq.(c-3),/,+,^ 

A  natural  requirement  of  any  numerical  flux  function  is  that  it  should  reduce  to  the  exact  flux 
in  certain  simple  cases.  For  example,  if  the  solution  is  a  constant,  q-q,  then  the  flux  is  a  constant, 
f=fiq).  Therefore,/'  should  satisfy: 

= /{9)-  (C-4) 

Schemes  of  the  form  shown  in  Eq.  (C-2)  for  which  this  condition  holds  (and  for  which /  depends 
sufficiently  smoothly  on  its  arguments)  are  said  to  be  consistent.  In  particular,  the  central 
difference  example  mentioned  above  is  consistent  since JJ+i/i  {q",  q/*i)  =  l/2(Jf"+J5*i)  implies  that 
//+  i/2y  =l/2(f  ("9 )  +/  (^))  =f('q).  Note  that  this  particular  notion  of  consistency  is  an  ex¬ 

ample  of  a  general  notion  that  may  be  more  familiar  to  the  reader.  Specifically,  a  finite-difference 
equation  is  said  to  be  a  consistent  approximation  to  a  differential  (or  integral)  equation,  if  the 
solution  to  the  latter  satisfies  the  former,  except  for  a  residual  that  vanishes  in  the  limit  of  ever 
increasing  grid  refinement. 

C-3.0  TVD  PROPERTY 

Although  the  central  difference  example  given  above  easily  demonstrates  consistency,  it 
does  not  lead  to  a  stable  scheme.  In  fact,  it  can  lead  to  errors  that  increase  without  bound  as  the 
number  of  time  steps  increases.  On  the  other  hand,  there  are  schemes  that  give  bounded  solutions, 
but  with  spurious  oscillations  near  discontinuities.  When  such  oscillations  develop  in  the  grid 
function,  { 9/ } ,  its  total  variation, 

} 

increases  with  respect  to  n.  Observe  that  this  is  an  appropriate  measure  of  the  spatial  variation  in  q", 
since  TV({q.''})  approximates  the  spatial  integral  of  I  d//dx  I.  Thus,  a  particular  notion  of  stability 
based  on  the  total  variation  is  now  introduced.  Specifically,  a  method  is  said  to  be  total  variation 
diminishing  (TVD)  if 

TV({9;+n)<TV({?;}). 
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The  TVD  property  is  geometrically  straightforward;  nevertheless,  from  a  practical  perspective,  it  is 
important  to  have  criteria  that  can  be  verified  easily  to  show  that  a  given  scheme  is  TVD.  For  this, 
let  the  terms  of  the  scheme  be  rearranged  in  the  form: 


At 


+  Lj 


Ax 


+  l: 


Ax 


Ax 


+  Rj 


Ax 


(C-5) 


where 


A<  =  and  A,qr;  = 


To  make  this  process  more  concrete,  consider  an  upwind  example  for  which  the  numerical  flux 
function  is  given  by: 


f'M  = 


3  +  1 


+  /;•) 


-  2I 


j+ 


il(^“+i  -  where 


/T-ri  -  /" 


(C-6) 


Observe  that  the  consistency  criterion  in  Eq.  {C-4)  is  satisfied  by  this  example.  Also,  a  few  calcu¬ 
lations  show  that: 


where 


«’*  ,  -  |a\ , 
.(+5  '  y+f 


A.9;  +  i  +  !«;_!  I  v,,;* 


Sr/;  ^ and  = 


The  upwinding  character  in  this  formula  is  seen  by  noting  that  backward  or  forward  differencing  is 
used  on  g  depending  on  the  sign  of  a,  i.e.,  the  wind  direction.  When  this  expression  is  inserted  into 
the  conservation  form  in  Eq.  (C-2),  the  following  result  is  obtained; 


Atg7 

At 


A,qp_j 

Ax 


Compare  this  with  Eq.  (C-5)  to  find  that  for  the  upwind  scheme: 


0.  K*  =  i 


;+* 


—  a 


J+2 


and 


Thus,  as  this  example  illustrates,  a  given  scheme  can  always  be  written  in  the  form  shown  in  Eq. 
(C-5).  Furthermore,  provided  At/Ax  is  sufficiently  small,  the  upwind  scheme  is  seen  to  be  TVD 
according  to  the  following  general  result  due  to  Harten,  (Ref  63),  and  Jameson  and  Lax,  (Ref.  64). 
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A  fmxtc-d%ffeTE.n(x  scheme  in  the  form  of  Eq.  (C-5)  is  TVD  if  <  0, 
ij  >  0,  Rf  >  0,  RJ  <  0  and  -  R-)At/Ax  <  1. 

Note  that  the  conditions  stated  here  are  sufficient  for  a  scheme  to  be  TVD.  They  are  not  necessaiy . 
Moreover,  a  given  scheme  is  not  represented  uniquely  in  the  form  shown  in  Eq.  (C-5).  Therefore, 
it  may  be  necessary  to  experiment  with  different  ways  of  defining  the  coefficients  before  the  above 
conditions  can  be  verified. 

C-4.0  CONVERGENCE 

The  forgoing  concepts  are  now  summarized  precisely  in  the  following  convergence  result. 

Suppose  a  grid  function  is  generated  by  a  TVD  finite-difference  scheme  'which  is 
consistent  with  the  scalar  wave  equation  in  Eq.  { C-l).  Also,  assume  that  the  scheme 
can  be  written  in  conservation  form.  Then,  the  grid  function  approximates  some 
wetdc  solution  to  Eq.  (C-l)  (in  an  integral  sense)  with  arbitrarily  small  error, 
provided  At  and  Ax  are  sufficiently  smedL 

For  example,  the  basic  upwind  scheme  defined  by  Eq.  (C-6)  satisfies  all  stated  conditions.  There¬ 
fore,  it  can  be  used  to  compute  an  arbitrarily  accurate  approximation  to  a  weak  solution.  However, 
this  scheme  is  only  first-order  accurate  and  is  extremely  dissipative.  In  other  words,  it  does  not  offer 
high  resolution  of  a  weak  solution,  without  a  highly  refined  grid.  Therefore,  the  discussion  now 
turns  to  the  development  of  high-resolution  schemes  that  satisfy  the  conditions  of  the  convergence 
result  given  above. 

C-5.0  HIGH-RESOLUTION  SCHEMES 

The  construction  of  a  high-resolution  method  begins  with  the  observation  that  the  order  of 
the  method  should  switch  according  to  the  smoothness  of  the  solution.  Specifically,  greater 
accuracy  can  be  achieved  in  regions  where  the  solution  is  smooth,  by  using  a  higher-order  method. 
On  the  other  hand,  such  a  method  causes  oscillations  around  discontinuities  such  as  shocks.  Thus, 
in  these  regions,  it  is  necessary  to  switch  to  a  lower-order  method. 
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C-6.0  FLUX-LIMITED  APPROACH 


This  switching  can  be  accomplished  according  to  a  flux-limited  approach  such  as  the 
following.  Let/'  be  a  numerical  flux  function  such  that  approximates/  to  an  order  cor¬ 

responding  to  the  smoothness  of  q.  In  other  words,  8 ^*/Ax  switches  from  a  high-order 
approximation  to/  where  q  is  smooth,  to  a  low-order  approximation  to/  where  q  is  nonsmooth. 
Specifically,/*  can  have  the  form 


= 


(C-7) 


where  8^  ‘*/A  x  and  hj  ^/Ax  are  high-  and  low-order  approximations  to/,  respectively.  Also,  ^  is 
a  nonlinear  flux  limiter.  It  is  so-named  because  the  use  of  high-order  flux  terms  is  limited  as  O 
switches  from  one  to  zero  whenever  a  lack  of  smoothness  is  detected  in  q. 


For  example,/”  might  correspond  to  central  differencing  as  explained  after  Eq.  (C-3).  This 
can  be  implemented  by  setting 


(C-8) 


Also,  /  *'  might  correspond  to  upwinding  as  explained  after  Eq.  (C-6).  This  can  be  achieved  with 


(C-9) 


For  the  flux  limiter,  <I>,  to  perform  the  switching  from/”  to  it  must  be  constructed  to 
detect  a  lack  of  smoothness  in  the  solution,  q.  For  example,  this  can  be  accomplished  by  monitoring 
the  ratio  of  differences. 


-  % 


=  0, 


^71  _ 

<lj 


Specifically,  r"  is  negative  when  the  solution  is  oscillating.  This  should  trigger  the  use  of  only  low- 
order  flux  terms.  On  the  other  hand,  when  r  "  is  close  to  one,  q  is  smooth  and  the  high-order  flux 
terms  should  be  switched  on.  A  typical  limiter  can  be  implemented  by  deflning: 


(C-10) 


where 


(p{r)  =  iiiinmod(  l,r) 


f  inin(  1 ,  r)  r  >  0 
[  0  r  <  0 


and  •0(r)  =  <^(i).  (C-Il) 
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Note  that  becomes  zero  when  r“<  0  or  j  <  0;  it  approaches  one  when  r"  and  |  are 
close  to  one.  This  gives  a  continuous  switching  required  in  Eq.  (C-7).  Also,  4>  can  be  shown  to 
satisfy  conditions  which  guarantee  that  the  scheme  is  TVD  (Ref.  65). 

C-7.0  SLOPE-LIMITED  APPROACH 

Alternatively,  the  switching  can  be  accomplished  according  to  a  slope-limited  approach  as 
follows.  As  before,  the  numerical  flux  function  is  constructed  so  that  5 /Ax  approximates/^  to  an 
order  corresponding  to  the  smoothness  of  q.  However,  in  the  present  approach,  the  switching  occurs 
among  q  values  instead  of  among  fluxes.  This  is  performed  within  the  framework  of  a  discretization 
procedure  which  is  originally  doe  to  Gudonov  (Ref.  66),  which  is  described  in  detail  below. 
Basically,  the  procedure  consists  of  the  following.  At  a  given  time  level,  the  variations  in  the 
solution  are  approximated  by  jumps  at  cell  interfaces.  In  fact,  the  accuracy  of  this  approximation 
determines  the  order  of  the  method.  Next,  the  waves  resulting  from  these  jumps  are  propagated 
forward  in  time.  When  this  information  reaches  the  next  time  level,  it  is  averaged  to  complete  the 
time  step. 

To  implement  this  procedure,  the  solution  is  initially  approximated  by  a  constant  on  opposite 
sides  of  a  cell  interface.  Let  the  value  on  the  right  of  x^iyj  be  denoted  by  q^+1/2,  and  the  value  on  the 
left  by  example,  they  can  be  determined  by  zeroth-order  extrapolation  according  to: 

7^='/“+!  =  (C-12) 

Since  these  values  are  generally  not  equal,  this  piecewise  constant  approximation  creates  a  jump, 
^*■1(2  ■  *be  cell  interface,  x^iq.  Next,  this  information  is  propagated  forward  in  time  by 

assigning  it  as  the  initial  state  for  the  linearized  problem. 


where 


XH  _ 


This  Riemann  problem  can  be  solved  explicitly  for  every  interface  Then  the  solutions  are 
superimposed  so  that  the  grid  value,  ^"^i,  is  determined  according  to  an  average  of  values  at  time 
level  r*'. 
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The  details  are  not  provided,  but  it  can  be  shown  that  this  scheme  has  a  numerical  flux  function  with 
the  form: 


(C-13) 


Note  that  this  is  identical  to  the  low-order  upwind  flux  in  Eq.  <  C-9),  under  the  condition  of  the 
zeroth-order  extrapolation  shown  in  Eq.  (  C-12).  This  gives  the  low-order  method  to  be  used  where 
the  solution  is  nonsmooth.  A  higher-order  method  is  obtained  by  improving  the  accuracy  of  the 
extrapolation.  However,  the  flrst-order  extrapolations, 


lead  to  unwanted  oscillations  in  nonsmooth  regions.  Thus,  to  limit  the  slope  variations  appearing  in 
these  linear  approximations,  slope  limiters,  ()>  and  >|/,  are  introduced  to  give: 

=  9^+1  -  ftj+l  =  </”  +  )V*9;- -  (C-14) 


Again,  (])  and  \|/  are  defined  by  Eq.  (C- 11),  but  their  interpretation  here  is  different.  Before,  they  were 
used  to  limit  the  use  of  high-order  flux  terms.  Here,  they  limit  the  use  of  large  slopes  in  a  linear 
approximation  of  the  solution. 


The  scheme  will  now  be  summarized.  When  a  lack  of  smoothness  is  detected  in  the  solution, 
the  interface  values  are  determined  by  zeroth-order  extrapolation.  This  leads  to  a  low-order  upwind 
numerical  flux  function.  On  the  other  hand,  when  the  solution  is  smooth,  the  interface  values  are 
determined  by  first-order,  linear  extrapolation.  The  result  is  a  high-order  upwind  numerical  flux 
function.  The  switching  is  performed  according  to  the  formulas  in  Eqs.  (C-13)  and  (C-14).  Also, 
the  slope  limiters  can  be  shown  to  satisfy  conditions  which  guarantee  that  the  scheme  is  TVD.  This 
slope-limited  method  is  an  example  of  a  MUSCL  scheme  (Ref.  67). 


C-8  VECTOR  FORM 


Now  consider  the  generalization  of  the  above  schemes  for  non  scalar  problems  such  as  the  one- 
dimensional  Euler  equations, 

Ql  +  fr  =  0. 
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Focus  on  generalizing  Eqs.  (C-9)  and  (C-13),  in  particular,  since  it  is  not  clear  how  to  form  a 
counterpart  to  aj^  ta  or  1/2.  This  is  accomplished  in  the  following  first-order  upwind  Roe 
scheme: 

^  +  ^  =  0 
At  Ax 

where 

/^v.  +  -  Qy). 

For  this,  the  absolute  value  of  a  matrix  is  defined  in  terms  of  the  absolute  values  of  its  eigenvalues. 
Specifically,  if  i4  =  S'^AS,  where  A  =  diag{X..},  then 

\A\  =  .S'" ’lAI.?,  where  |A|  =  tliag{|A,|}. 

The  most  conspicuous  property  defining  above  is  that 

Such  a  matrix  can  be  constructed  by  a  special  averaging  procedure.  TTiis  was  done  by  Roe  (Ref.  68) 
for  the  case  of  a  perfect  gas,  and  extended  to  nonequilibrium  flows  by  Liu  and  Vinokur  (Ref.  69). 

The  flux-limited  approach  introduced  above  can  be  implemented  by  setting 

=  ^{^Ai  +  /^;*)  -  -  Qn 

where  $"+ 1/2  is  a  diagonal  matrix  whose  ith  entry  is  a  flux  limiter  as  shown  in  Eq.  (C-10),  with  q" 
replaced  with  the  ith  component  of  SD]. 

Finally,  the  slope-limited  MUSCL  approach  introduced  above  can  be  implemented  by 

setting 

where 

=  F(Qf^p  - 
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and  the  components  of  Qf^,  and  are  defined  in  terms  of  slope  limiters  as  shown  in  Eq. 

(C-14). 

Note  the  complexity  of  defining  such  upwind  schemes.  In  addition,  extending  upwinding 
techniques  to  higher  dimensions  is  not  straightforward.  In  fact,  attempts  to  generalize  these  schemes 
to  two  or  three  dimensions  have  resulted  in  the  production  of  zigzag  shocks  unless  the  grid  is 
aligned  with  the  shock.  Also,  the  eigenvalues  and  eigenvectors  of  A  would  have  to  be  rederived  and 
coded  whenever  a  new  thermo-chemical  model  was  considered.  Therefore,  as  an  alternative  to  the 
matrix-based  dissipation  models  presented  above,  Jameson's  flux-limited  scalar  dissipation  was 
selected  for  incorporation  into  NEDANA  to  satisfy  the  objectives  of  this  developmental  effort. 
Jameson’s  model  is  discussed  below. 

C<9.0  JAMESON’S  FLUX-LIMITED  DISSIPATION  MODEL 

The  high  resolution  schemes  based  on  upwinding  presented  in  the  previous  subsection  have 
developed  in  parallel  with  those  based  on  artificial  dissipation  models.  The  concept  of  artificial 
dissipation  originated  with  von  Neumann  and  Richtmyer  (Ref.  70),  who  were  attempting  to 
simulate  computationally  the  propagation  of  shock  waves  in  inviscid  fluid  flow  without  generating 
mesh  scale  numerical  oscillations,  von  Neumann  proposed  that  to  suppress  these  oscillations,  the 
difference  equations  could  be  augmented  with  terms  reminiscent  of  the  viscosity  terms  in  the 
Navier-Stokes  equations.  However,  the  proposed  artificial  terms  were  purely  numerical  and  did  not 
correspond  to  any  physical  dissipative  mechanism.  This  concept  has  evolved  so  that  modem 
implementations  operate  adaptively.  In  the  more  recently  developed  schemes,  smoothing  terms  are 
made  to  dominate  in  the  vicinity  of  discontinuities  when  a  lack  of  smoothness  is  detected  in  the 
solution.  In  particular,  Jameson’s  flux-limited  artificial  dissipation  model  has  the  simplicity  of  the 
artificial  dissipation  approach  while  satisfying  the  TVD  criteria. 

Subsection  6.0  gives  a  general  description  of  the  flux-limited  approach  to  achieving  the  TVD 
property.  Here,  the  particular  method  used  in  NEDANA  is  explained.  This  flux-limited  dissipation 
approach  was  Urst  studied  by  Jameson  (Ref  38),  and  modified  by  Yoon  and  Kwak  (Ref.  39),  and 
later  by  Deese  and  Agarwal  (Ref  71).  It  can  be  described  easily  for  the  scalar  Eq.  (C-1),  and  its 
finite-difference  approximation  in  Eq.  (C-2).  Here,  the  numerical  flux  function  is  defined  by 

/;+i  =  +/r’)  +  (C-15) 

where  <|)"+  /  =  <|>(i'"+*i*)  and  ‘)  ate  defined  in  Eq.  (C-11).  Also,  e^,y2  is  an  adjust- 

ed)le  parameter  which  can  be  used  to  control  the  amount  of  artificial  dissipation. 
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Recall  from  subsection  C-6.0  that  if  oscillations  are  detected  in  the  solution,  the  numerical 
flux  function  shown  in  Eq.  (C-7)  becomes  the  flux  shown  in  Eq.  (C-9).  Here,  oscillations  cause  4|) 
and y  to  vanish  and/j^.|/2  becomes 


This  resembles  the  low-otder  form  shown  in  Eq.  (C-9),  and  it  is  useful  to  view  the  second  term  as 
artificial  dissipation.  Next,  recall  that  if  the  solution  is  smooth,  the  numerical  flux  function  in  Eq. 
(C-7)  becomes  the  flux  shown  in  Eq.  (C-8).  Here,  when  q  is  smooth,  (|)  and  \|t  are  close  to  one,  and 
/j+i/2  becomes 


,n+l 


This  resembles  the  high-oider  form  shown  in  Eq.  (C-8)  since  the  second  term  is  negligible  in  terms 
of  truncation  enx»r;  thus,/**  here  and  in  Eq.  (C-8)  leads  to  an  approximation  of^^of  the  same  order. 
The  second  term  here  is  used  because  it  has  been  found  to  inhibit  the  odd/even  decoupling  of  grid 
values  that  can  result  from  the  use  of  a  central  difference  scheme.  Finally,  note  that  the  method 
defined  by  Eqs.  (C-2)  and  (C-1 5)  can  be  shown  to  be  TVD,  provided  Ej+ia  S  1.  The  numerical  flux 
function  for  Jameson’s  scheme  takes  the  following  form  for  systems  of  equations: 


(C-16) 

where 

cj  =  ^  A”+i]  ,  <‘+1  =  I'"'-*  +  «4inax(t/7,i^;Vi)]>  = 

P?+,-2)>7  +  p^, 

P^i+2pJ  +  P"-i 

Here,  'F  is  defined  for  vectors  P  and  Q  to  have  components: 

=  ^[sign(P.)-|-sign(Q,)]niin(|/’i|,|(3,|). 

Note  that  this  is  a  convenient  way  of  expressing  the  limiters  since  for  the  scalar  case, 

and  A^gJ+i). 

The  other  quantities  in  Eq.  (C-1 6)  include  and  which  are  adjustable  dissipation  parameters, 
and  X  which  is  the  spectral  radius  of  the  flux  Jacobian  dF/dQ.  An  important  advantage  of  the  present 
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method  is  that  only  the  spectral  radius  of  the  flux  Jacobian  need  be  rederived  and  implemented 
when  new  thermo^hemical  models  or  new  physical  phenomena  such  as  magnetohydrodynamics 
are  considered  for  inclusion  in  the  code. 

Consider  the  terms  that  comprise  First,  E^i/i  depends  on  the  pressure  sensors,  v"  and 
V  .  These  sensors  allow  the  dissipation  to  be  adjusted  adaptively.  Specifically,  they  are  large  in 
the  presence  of  high-pressure  gradients  and  negligible  in  regions  of  a  smoothly  varying  pressure. 
The  influence  of  these  variations  on  is  attenuated  or  amplified  according  to  the  value  of  K^. 
Also,  K2  is  a  constant  chosen  large  enough  to  suppress  small-scale  background  oscillations. 

This  implementation  of  the  flux-limited  dissipation  scheme  is  a  componentwise  application 
of  the  scalar  constructions  to  the  vector  equations.  However,  such  a  generalization  is  not  uniquely 
determined.  For  example,  the  quantity,  appearing  in  the  scalar  development,  has  a  natu¬ 

ral  matrix  counterpart  here.  Yet,  for  ease  of  implementation,  it  is  replaced  by  the  scalar  coefficient, 
c^.  Therefore,  this  is  called  a  scalar  dissipation  model.  The  form  of  the  numerical  flux  function  used 
for  the  quasi-one-dimensional  set  of  conservation  equations  may  be  found  in  Ref.  17. 
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APPENDIX  D 

DEFINITION  OF  TIME  STEP 


The  computational  time  step  at  the  cell  volume  J,  K,  L  is  set  to  be  the  minimum  of  the  time 
steps  in  the  individual  computational  directions 


where 


and  CFL  is  the  Courant-Friedrichs-Lewy  number.  The  spectral  radii  ate  defined  as 


where 


AfJ.ft-,/,  = 

a  •  fQ 

J.K.L  + 

\j,K.L 

Af,  J,K.L  = 

ii  ■  (tvj 

J.KJ.  + 

>/<iA  fij 

Ac  J,h\L  = 

u  - 

JJ<.L  + 

2 
1 


^  ['^c^A.r+1  + 

The  frozen  speed  of  soimd,  a,  is  calculated  from 

2  P 
a  =  V- 

where  the  frozen  ratio  of  specific  heats,  yj-,  is  expressed  as 

7/  =  (l+/i), 


(D-1) 


=  CFL 

V 

J.A'.t 

=  CFL 

’v' 

X 

1 

J,K,L 

=  CFL 

v' 

J,K.L 

(D-2) 

(D-3) 


(D-4) 


(D-5) 


(D-6) 


where  ^  =  dp/dE.  The  evaluation  of  p  is  dependent  on  the  type  of  nonequilibrium  model  employed. 
The  current  version  of  the  NEDANA  flow  solver  is  limited  to  NEQPAK  thermodynamic  models 
one  or  two.  Therefore, 


3 


_  I  p/T,  model  =  1 

P3(^'t,lr)  ip  -  Pe)/T,  lliodel  =  2 


(D-7) 
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APPENDIX  E 

QUASI-ONE-DIMENSIONAL  NOZZLE  FLOW  EQUATIONS 

The  purpose  of  this  appendix  is  to  derive  the  partial  differential  equations  that  model  inviscid 
quasi-one-dimensional  nozzle  flow.  The  derivation  proceeds  by  integrating  the  differential  form  of 
the  equations  over  a  nozzle  cross  section  of  vanishingly  small  width.  Then,  the  integral  is 
transformed  in  steps  by  applying  certain  assumptions.  For  example,  it  is  assumed  that  there  are  no 
azimuthal  variations  in  the  flow.  Also,  for  simplicity,  nozzle  cross  sections  are  assumed  to  be 
circular.  Finally,  the  required  result  is  obtained  in  the  limit  of  decreasing  cross-sectional  width.  See 
Ref.  72  for  more  information. 

To  facilitate  the  integration  of  the  differential  form  of  the  equations  over  a  nozzle  cross 
section,  the  full  three-dimensional  equation  set  is  first  expressed  in  the  curvilinear  coordinates, 

^  =  ^(Tr.j/.z),  7j=  i]{x,y,z),  C  = 

For  convenience,  the  notations 

{x,y,z)  =  (a;i.  =  x  (C,V,C)  =  =  ( 

{u,v,w)  =  (i>,.  =  iJ  {F,G,H)  = 


are  used.  Here,  F,  G,|md  H  are  the  Cartesian  flux  vectors  defined  in  Eq.  (46).  The  Jacobian  of  the 
transformation  Jc  |  is  written  as 


J  =  clet 


Following  Ref.  37,  the  differential  form  of  the  conservation  equations  is  written  in  the  curvilinear 
coordinates  as: 


<)Q  dF  dc 
dt  di]  dC 


(E.1) 
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where  for  flows  in  thermo-chemicai  nonequiiibrium 


Pi 

Pit' 

Pn.s 

PntV 

Ei> 

II 

EvU 

pu 

puU  +  pdi/dx 

pn 

pvU  +pd^/dy 

pw 

pvfU  +  pd^jdz 

E 

PiV 

P^W 

Wl 

PnsV 

Pn,W 

EvV 

.  H  =  T^ 

EvW 

u)y 

puV  +  pdrijdi 

puW  +  pdC,jdx 

0 

pvV  +  pdr)/dy 

pvW  +  pdC,fdy 

0 

pwV  +pdti/ds 

pwW  +  pdQfdz 

0 

{E-\-p)V 

[E  +  p)W 

0 

{/  =  ■  jj,  V'  =  •  iT.  W'’  =  VC  •  iT. 


Here,  ((/,  V,  W)  are  the  sO'Called  contiavariant  components  of  velocity.  For  convenience,  the 
following  notation  is  used; 


Then, 


([/,  V,W)  =  ('V', ,  Vi.  v:-})  =  V  and  {f,  (i,  J? )  =  (fi,  /’j,  /a). 


;=!  j  =  l 


1  <  i  <  3. 


(E-2) 


Now,  a  particular  nozzle  coordinate  system  is  introduced.  Let  the  Cartesian  coordinates  x,  y, 
and  z  be  situated  so  that  the  x-axis  is  aligned  with  the  central  axis  of  the  nozzle.  Then,  let  the 
curvilinear  coordinates  be  defined  so  that  ^  varies  only  along  the  length  of  the  nozzle,  i.e.,  it 
depends  only  on  x.  Also,  within  a  transverse  planar  section,  let  ^  be  a  radial  coordinate  and  an 
azimuthal  angle.  In  particular,  note  that  the  nozzle  surface  is  not  a  ^ = constant  surface.  Specifically, 

C  =  C(j‘).  V  =  C  =  + 

y  —  (^  COST].  z  =  — C  sin  tj 
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and 


J  —  det 


r  ^ 

^  1 

dx 

dy 

dz 

£2 

dx 

dy 

dz 

K 

K 

££ 

.  dx 

dy 

dz  . 

=  det 


0 

0 


0  0 
-C”’sinT/  -^~*cosf; 
cost;  —  sini; 


-  e  /--i 


(E-3) 


The  quasi-one-dimensional  equations  are  derived  by  integrating  the  three-dimensional  equations 
over  a  cross-sectional  volume,  V,  defined  by; 


V  =  ^0  <  4  <  fo  +  At,  0  <  7/  <  277,  0  <  C  <  r(t)}  , 


where  ^  =  r(4)  defines  the  nozzle  surface.  This  leads  to 

t'iT  niO  r .  .  .  .  , , 

/  /  /  0(  +  +  fiij  +  d(^d^di]  =  0. 

Jo  Jo  ^  ^ 

Now  using  Eqs.  (E-2)  and  (E-3)  and  the  lack  of  azimuthal  variation  in  each  Fj, 


d 

dv 


j=l  ^^3 


(&4) 


Next,  for  fixed  1),  the  following  is  an  integration  over  an  area,  say^l,  in  a  meridian  plane  of  the 
nozzle.  Green’s  theorem  in  the  plane  gives: 


-H,F)fdT 


(E-6) 


where  d,4denotes  the  boundary  of  the  area,>l,  ^  denotes  a  unit  vector  tangent  to  dA  with  ,4  to  the 
left  oft,  and  dx  denotes  an  infinitesimal  line  element.  On  the  edges; 

^  =  ^0  +  Af,  (0.  1).  dr  =  dC,  ^  =  ^0,  T  =  (0,-1),  dr  =  -dC: 

C  =  r(0,  f  =  .  dr=  -Jl  +  (re)2df;  C  =  0,  f  =  (1,0),  dr  =  df . 

yi  +  (t-£)^ 
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Thus, 


/“  -  -  I  /■«o+A{  „  ,  , 

+  /,  •{0,1K+/ 


-I 


’■(^o+Af)  . 

F 


/■’■(to)  -  /40+A(  .  I 

-  F  dc.-  h\  d^. 

Jo  i=(.o  J^o  'C=« 


«=6)+Af 


^o+A(  . 


(E-7) 


The  last  term  vanishes  since  by  Eq.  (E-3), 


h\ 

!<=0  ^  Ox, 

7—1  ■' 


i-V  F 

'dx, 


oc 


c=0 


=  0. 


(E-8) 


C=o 


Next,  for  the  second  term  on  the  right  side  of  Eq.  (E-7),  it  will  be  shown  that 


P^(-T^V  +  W) 

0 

P..A-nV^w) 

0 

Ev[-r^U  +  W) 

0 

pn{-r^l!  -I-  M^)  - 

-rpr^ 

pi'{-r^U  -1- TV) -I- p cost; 

'■fJ^COST/ 

pw{  —Tf^U  -|-  W)  -  p sin  1} 

-r^j’psiniy 

{E-Vp){-T^U  ■\-W) 

C=r 

0 

(E-9) 


The  first  equality  here  follows  from  Eq.  (E-3)  and  the  definition  of  P"  and  ^ .  Also,  the  fact  that 
(-r^(/  +  IV)  is  zero  at  the  nozzle  surface,  where  ^  =  r(^),  can  be  determined  by  ^plying  a  tangent 
flow  boundary  condition.  For  this,  note  that  the  nozzle  surface  is  represented  by/=  ^  -  r(4)  =  0. 
Therefore,  the  outwardly  directed  unit  normal  is: 


_  (->’£^3,,COST/,-sin7?) 
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Also,  the  contravariant  velocity  components  ane  given  by 


U 

V 

w 


=  uf* 

ii-Vi)  =  — ^“’(usin  7/ +  lii cost;) 

v-V(  =  u cos i;— nisi n  t;. 


Thus,  the  vanishing  of  the  normal  component  of  velocity  at  the  nozzle  surface  means  that; 

n  -  +  u cost;  -  i/i sin  t;  -rf[!-^-W 

0  =  u  ■  n  =  - i =  -j-  >  =  0. 

Now  Eq.  (E-9)  follows  from  this  formula.  Finally,  for  the  first  and  third  terms  in  Eq.  {E-7),  note  that 
by  Eq.  (E-3)  and  the  fact  that  ^  is  a  fiinction  only  of  x. 


Combining  Eqs.  (E-4 )  -  (E-10)  and  making  use  of  Eq.  (E*3  )  gives; 

/•2ir  f  fio+M  /’•(^)  r  frlio+^i) 

0=/  u  /  ^[Qi-a]dc  d^+ 

Jo  sx  J  JO  ' 


(E-10) 


d( 


r(ta)  [ 


<0 


0 
0 

-rpr^ 

’’fr  Vcost; 
-r^~'psiii  t; 

0 

According  to  the  azimuthal  independence  of  Q,  F  and  Q, 

r  frU)  /■  1  /’■(io  +  Af) 

0  =  2ir/  /  -^[Qt-^UC  d^  +  2n 

Jio  Jo  sx  JO 


jdi;. 


rf 


0 

0 

-rpr^ 

0 

0 

0 


di. 
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Observe  that  certain  components  vanish  in  the  last  term  since  cos  sin  s  0.  Next, 

this  equation  is  transformed  by  applying  the  following  mean-value  theorem  (Ref.  73).  Suppose  that 
/and  g  are  integrable  on  a  set,  S,  where  g  ^  0.  Iben,  there  is  an  average  value,/  between  min  ^/and 
max  /  such  that 

ImsKW  =  f 


This  result  can  be  applied  to  the  above  equation  componentwise.  For  this,  define  the  area  function, 


=  27r  /  so  that 

Jo 


A(^  =  27rrr^. 


Letting  A(^)  play  the  role  of  J  in  the  above  mean- value  theorem  gives; 


“ = C  [(^«),  -  (u“)] + II 


0  1 

0 

0 

■A^p 

0 

0 

0 


where  Q ,  F  and  Oi  here  represent  average  values  over  a  given  cross  section.  Dividing  by  A^,  and 
taking  the  limit  as  ^  0  leads  to: 


r  0 


{C^aq)^  +  {af)^=  (cMft)-h 


0 

0 

A^p 

0 

0 

0 
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Since  the  equations  involving  pv  and  pw  are  decoupled  from  the  others,  the  nozzle  flow  equations 
are  taken  as 


\ 

/ 

P\<x 

AJ  '  wi 

■AJ-^ 

Pns 

AJ-' 

PnsU^x 

Ev 

EvU^x 

AJ~^  liv 

pu 

(/>“'■*  +  P)^x 

AiP 

\ 

E 

1 

\ 

.  (E  +  p)u^j.  . 

) 

i 

0 

where  7  =  ^  is  the  Jacobian  of  the  transformation  x  Observe  that  by  using  the  fact,  y  1 1 
certain  terms  above  can  be  cancelled.  However,  this  form  is  retained  to  parallel  the  three- 
dimensional  case  in  Eq.  (E-1).  The  finite-volume  form  of  the  equations  is  easily  obtained  by  taking 
J'^A  as  the  volume  of  the  corresponding  cell  (Ref.  37).  Finally,  note  that  setting  A  =  1  here  gives 
the  equation  set  for  a  shock-tube  problem. 
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APPENDIX  F 

THE  NONEQUILIBRIUM  SOURCE  JACOBIAN 

F-1.0  DEFINITION  OF  SOURCE  JACOBIAN 

The  numerical  scheme  is  based  on  the  linearization  of  the  nonequilibrium  flux  vector.  The 
vector  of  conserved  variables,  Q,  and  the  nonequilibrium  source  vector,  O,  for  the  NEDANA  flow 
solver  are  defined  as  follows 


,  Q  = 


The  chemical  source  terms,  to,,  are  the  production  of  species  s  in  units  ig/(m^s).  The  nonequilibrium 
energy  sources,  (0^„  are  the  production  of  energy  i  in  units  The  nonequilibrium  source  vector 

is  linearized  as  follows: 

=  Sr  +  Z”  C7(Af2),  (p.2) 

where  Z"  is  the  Jacobian  of  11"  with  respect  to  |2">  and  5|2"  =  -  Gt-  '1^6  matrix,  Z,  has  the  form 


3pi 

0 

0 

0 

0 


9pn-  OEn  I 
0  0 


_'^wi  d^f 

BEn 


gu*tiji  duhta 

3J^uIy  9E 

du/Ki  dufNt  SufWl 

^^ffnr  ^(^*0  c?(ptr)  SJpun  8E 


Swrifit 

9Effnt  3(pu)  d\pv)  d{pw)  9E 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 
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The  scope  of  the  cunent  flow  solver  is  to  use  full  finite-rate  chemistiy  with  a  two- 
temperature  thermal  model.  In  this  model,  the  translational  and  rotational  energies  are  assumed  to 
be  in  equilibrium  with  one  another  at  the  temperature,  T.  The  vibrational  energy  is  allowed  to  be 
characterized  by  a  separate  temperature  Ty.  In  the  case  of  ionization  or  electronic  excitation,  the 
electron  temperature,  T„  is  assumed  to  be  in  equilibrium  with  the  vibrational  temperature  at  the 
temperature,  The  nonequilibrium  energy  modes  of  all  species  are  characterized  by  a  single 
temperature.  Therefore,  the  current  scheme  has  only  one  nonequilibrium  energy,  Ey.  This  energy 
contains  the  vibrational,  electron,  and  electronic  energies  of  all  species.  The  NEDANA  flow  solver 
employs  NEQPAK  to  provide  the  nonequilibrium  source  terms  and  their  derivatives.  The 
adulation  of  NEQPAK  to  the  development  of  the  sources  and  their  derivatives  will  now  be 
discussed. 

F-2.0  DERIVATION  OF  CHEMICAL  JACOBIANS 
Define: 


II 

(F-4) 

(F-5) 

713 

p  =  ps, 

5=1 

(F-6) 

ns 

.-w  =  E  - 

7^1  T' 

(F-7) 

II 

{F-8) 

U9 

=  E 

3=1 

(F-9) 

ns 

(F-10) 

a=l 


where  Y,  is  the  concentration  of  species  s  in  kgmoleln?\  y  is  the  concentration  of  the  mixture  in 
kgmolein?',  is  the  density  of  species  s  in  kgln^\  p  is  the  mixture  of  the  density  in  kgln?\M.^  is  the 
molecular  weight  of  species  s  in  kg/kgmole\M  is  the  molecular  weight  of  the  mixture  in  kg/kmole, 
and  y,  is  the  mass  fraction  of  species  s.  Cy„  is  the  specific  heat  at  constant  volume  dependent  on  T. 
C,y  is  the  specific  heat  at  constant  volume  dependent  on  Ty.  Also,  define  a  general  temperature  array 
such  that 
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T„.  =  T„.(T.Tk),  ni  =  l.ntypc. 

For  example,  with  raype  =  2, 

Ti  =  T 
T2  =  T'~°T^ 


(F-11) 


(F-12) 


where  0<  a  <  1. 

NEQPAK  provides  {kmoletm^lsec)  which  is  a  function  of  y,  and  T„.  NEQPAK  also 
provides^  .^,and  ^  ,  where  the  vertical  bars  denote  the  partial  derivatives  are 

evaluated  holding  the  subscripted  quantttes  constant.  The  task  now  is  to  write  these  quantities  in 
terms  that  the  NEDANA  flow  solver  requires.  The  chemical  sources  become 


«>i‘  =  Mi  Wf 


(F-13) 


The  total  derivative  of  wj  becomes 

dii>,  =  -Vi ,  dti>, 

dij. 


-  vi  (  V 


ntvpe  ... 

+  y'  — - 

dT„, 

tll=1 

duj 

«l=l 


dyj 

It. 


(F-14)' 


The  total  derivatives  d^,,  dT,  and  dTv  must  now  be  expressed  in  terms  of  the  conservative  variables, 
Q.  First, 

fhj  =  (F-15) 


Now,  consider  dT.  The  total  energy  of  the  mixture  is 


£’  =  ^/)(  +  tiJ^)  +  £;  +  £k- 

The  internal  energy  dependent  on  T  is 


(F-16) 


£;  =  Pf/. 


(F-17) 
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where 


and 


e/  = 

1=1 


(F-18) 


(F-19) 


where  C|,_  ,r  and  ef  are  the  specific  heat  at  constant  volume  of  species  (  due  to  T,  and  the  energy 
of  formation  of  species  i,  respectively.  Then, 


(kj 


ns  ns 


1=1 


1=1 


where 


-  C' 

—  ^  II. ir 


dT 


(F-20) 


(F-2I) 


Thus, 


Ens  V 

1=1  '  v,tr 


(F-22) 


Now, 


dY, 


dp,  -  y,  dp 
P 


and 

ns 

f^p  = 

■=i 

Now,  writing  E  in  terms  of  Q, 


E  = 


{puf  +  {pvf  +  ipw^ 


and 


(F-23) 


(F-24) 


(F-25) 


dJS  =  ud{pti)  vd[pv) wd[pw)  pdej d{Ev) 
+[/■/  -  +  ij*  +  w'^)]dp. 


(F-26) 
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Using  the  fact  that 


r/  =  t  -  -  fcV/, 


(F-27) 


substitution  gives 

dE  =  +  w»d(/ow)  +  prfe/ +  d(£v) 

+  [r  -  (Ki'^  +  -  ry]dp. 


(F-28) 


Rearranging, 

dci  =  -|ri(£’)  —  ud{pu)  —  vd{px))  -  wd{pw)  -  d{Ev) 

-[f  -  ( ti^  +  —  ev]tip|.  (F-29) 


Now, 


dT  =  — <  d{E)  —  ‘tid{pu)  —  vd{pv)  —  wd(pv.i)  —  d{Ev) 

na 

-[r-  (7i’^  +  t;^  +  !/»*)  -  ev]dp  -  ^  p/.,[<ip,-  -  Yidp]^  (F-SO) 

>=i 

leading  to 

dT  =  ■— — |f/(£)  -  ud{pu)  —  vd{pxi)  -  v>d(pw)  -  d{Ev) 

P^  Utir  '* 
ns  I 

.=1  ^ 

Now,  to  express  d7V  in  terms  of  the  conservative  variables  Q.  The  vibrational/electronic  energy  is 


Ey  =  pfV'  (F-32) 


where 


and 


7i9 

•'v  = 

1=1 


(F-33) 


(F-34) 
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where  ^  the  specific  heat  at  constant  volume  of  species  i  due  to  Ty.  Then, 


drv  =  ^  dY^  €v,,  +  £  V;  dev,.- 


(F-35) 


t=i 


1=1 


Thus, 


But  also, 


Therefore, 


Now,  define 


_  dev  -  E:=i  dYey,,  _  dey  - dY, 
2^i=\  t 


d{ev)  = 


dEy  —  cydp 


u,V  [ 


d[Ey)  -  Pv.,  dp, 
1=1 


■1 


(F-36) 

CF-37) 

(F-38) 


dCv, 

W 

t?w, 


du, 

iri=l 

du^i 


dT,. 


-t,7 


m=l 


aT„ 


dT,.. 


dTy 


(F-39) 


Employing  Eqs.  (F-15),  (F-31),  and  (F-38),  Eq.  (F-14)  can  now  be  written  in  terms  of  derivatives 
of  the  conservative  variables. 


Mi  dtj, 

pC„,tr  dT 


+ 


I'l.Ty 


[ti(£')  -  ud{pu)  —  vd{pv)  — 


■  Af. 

Af, 

dd}. 

pC\,,y 

dTy 

'V.7’ 

P^\',lr 

w 

■7,T^' 

+y'(  — 


M, 
,  T  Mj 


+ 


Aft 

pC„,tr  [‘2 


-(«'=+  f/  +  W^)-  e;,j 


dui 

w 


The  chemical  source  Jacobians  can  now  be  obtained  by  taking  partial  derivatives  of  the  above 
equation. 

rtiMi. 

=  +- 

j  ,^U,ptpl,p2M,£^V 


diift 

dpj 


d{Ev) 

Mi  dH), 

•>,7V  pf^'v.v  dTy 

»1  by  takii 

Mi  dui 


7.T 


}dp} 


(F-40) 


Afj 


+ 


M, 

pC-V^tT 

Mi 

f^\iy 


-{li^  +  r;'^  «!■*)  -  e/., 


dCji 

W 


evj 


du>. 


y.Ty 


dTy 


->.T 


(F-41) 


124 


AEDC-TR-94'18 


duit  _  Mj  dHjj 

„.pu.iyv.pw.E  p(^'v,tT  (iT  -y.Tv 


M,  dCjt 


^pCv.V  dTvU^T 

(F-42) 

Mt  r)w, 

QT  'ijv 

{F-43) 

Ml  OCjy 

‘p.pn.pv;.Ev,E 

pCv.tT  dr 

(F-44) 

du)i 

My  OCiy 

d{pw) 

nPiTv 

(F-45) 

M% 

P,/niii<ii,i)'<i\Ev 

^  p(^'v,tT  dT 

(F-46) 

F.3.0  DERTVATION  OF  VIBRATIONAL/ELECTRONIC  JACOBIANS 

The  vibratLonal/electronic  source  term  for  the  two-temperature  model  where  Ty=  Ty=  T,  has 
the  form 

•"‘i  =  -/’l  l.  +  ^Vc  +  ~  }tf  V  -u.  (F-47) 

The  first  term,  is  the  Landau-Teller  relaxation  term  between  T  and  Ty.  NEQPAK  returns  this 
term  in  J/rn^/s,  as  well  as  the  following  derivatives:  j,  j,.  The 

third  term,  (0^,,  is  the  relaxation  term  between  T  and  T,.  NEQPAK  provides  this  term  in  7/(m^j), 
^  well  as  the  following  derivatives:  jv  ’  L  t'  t®*™"  Pe 

V  •  H ,  is  the  electron  pressure  gradient  term,  ihiis  term  is  not  provided  by  NEQPAK.  This  term 
is  treated  as  a  viscous  term  and  does  not  appear  in  the  Jacobian  formulation.  The  second  term, 
is  the  production  or  destruction  of  Eydue  to  chemical  reactions.  This  term  is  currently  not  provided 
by  NEQPAK.  Instead,  due  to  its  simple  form,  it  is  developed  in  the  NEDANA  flow  solver.  The 
non-preferential  form  is 

ns 

WVV  =  WfCv',,. 

■=1 
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or 

n« 

UVc  =  51  WiCV.t. 

t=l 


Using  Eq.  (F-39),  the  derivatives  are  easily  formed  to  yield 


clwt  - 


OuJVc 


y,i^].T.Tv 


u*  a* 

_  ^  OiiSi 
■ 


^V,t 


dT 


~i,Tv 


rti  sj“ 

dT 

t=l 


^Va 


du)vc 


dTv 


'\.T 


,  =  L 


1=1  L 


dui 


dTv 


1.T 


fiv.i  + 


(F-49) 


(F-50) 

(F-51) 

(F-52) 


The  derivatives  of  the  total  source  term  can  now  be  formed  by  combining  the  individual  terms  to 
yield 


diOv 


d-ij 


du>vi, 


dxi 


+ 


doJVf 


d^j 


+ 


durvc 


-„5ij,T.7V 


^7, 


^i,T.Tv 


(F-53) 


dujv 

'W 


f)uv„ 


f,Tv 


i)T 


+ 


dt^vt 


•yj'v 


dT 


nt.'/V 


+ 


duvc 


dT 


tt.Ti' 


dwv 

dwVr 

dwvc 

dwvu 

dTv 

-..T  "  ^Tv 

->.T  '^V 

.j  ^  dTv 

(F-54) 

(F-55) 


The  total  derivative  of  tOv  becomes 

dwv  = 


duJv 


Uljv 


dyj  + 


OTv 


dTv 


dT 


y.T 


(F-56} 


The  total  derivatives  d<f,  dT,  and  dT,,  were  defined  in  Eqs.  (F-15),  (F-31),  and  (F-38),  respectively. 
Using  these  definitions  of  the  total  derivatives,  Eq.  (F-S6)  becomes 
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dwy  = 


1  dufy 

pa,,tr  'W 


[f/(£^)  —  u(f(^u)  —  vd(/>ir)  — 


1 

1  du^v 

pC'.yV  OTv 

-I,T  pC'v.tT  dT 

Ljv. 

d(Eir) 


f  1 


T.^./.7'.7V 


+ 


1  n 


pCVtr  L2 


-(  „  +  „  +  )  _  c;_j 


duy 


OT 


evj  du>v 


',Tv  P^‘v,v  S'Tv 


(F 


57) 


The  vibrational/electronic  source  Jacobians  can  now  be  obtained  by  taking  partial 
derivatives  of  the  above  equation. 
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E 

f 

■* 

F 

AF 

F,G,H 


NOMENCLATURE 

Frozen  speed  of  sound,  mis 
Coefficient  in  Eq.  (13) 

Scalar  flux  Jacobian 
One-dimensional  nozzle  area,  n?' 

Matrix  flux  Jacobian 
Millikan  and  White  coefficient 
Arrhenius  coefficient,  Eq.  (30) 

Diagonal  matrix,  1/s 
Coefficient  in  Eq.  (20) 

Iteration  matrix 
Block  iteration  matrix 
Scalar  iteration  matrix 
Viscous  term  coefficient 
Arrhenius  coefficient,  Eq.  (30) 

Scalar  dissipation,  n?ls 
Park’s  correction 
Arrhenius  coefficient,  Eq.  (30) 

Specific  heat  at  constant  pressure,  Jlkg  K 
Specific  heat  at  constant  volume,  Jlkg  K 
Courant-Friedrichs-Lewy  number 
Coefficient  in  Eq.  (23) 

Effective  diffusion  coefficient,  mVs 
Binary  diffusion  coefficient,  mVs 
Specific  energy,  Jlkg 
Energy  per  unit  volume,  Jln^ 

Scalar  flux  function 
Numerical  flux  function 
Change  in  free  energy,  Jlkmole 
Cartesian  flux  vectors 
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F,  G,  HF 

Si 

G 

G 

h 

h 

ft? 

'i  'i  'i 

'jp  V 

I 

jXl 

J 

k 

kf 

L 

C 

M 

M 

M 

ne 

net 

nq 

nr 

ns 

P 

9 

9 

Q 


Computational  flux  vectors 

Degeneracy  of  /th  electronic  level 

Gibbs  flee  energy,  Jfkmole 

Rate  of  gain,  kmoletn^s 

Planck  constant,  Jts 

Enthalpy  per  unit  mass,  Jlkg 

Species  heat  of  formation,  Jlkg 

Cartesian  unit  vectors 

Identity  matrix 

Computational  indices 

Jacobian  of  transformation 

Boltzmann  constant,  JIK 

Forward  reaction  coefficient,  m^fkmole 

Reverse  reaction  coefflcient,  m^/kmole  or  m^lkmol^ 

Equilibrium  constant 

Characteristic  length,  m 

Rate  of  loss,  hnolein?  s 

Generic  molecule 

Mach  number 

Molecular  weight,  kg/kmol 

Avogadro’s  number 

Number  of  nonequilibrium  energies 

Number  of  electronic  energy  levels 

Number  of  conserved  variables 

Number  of  reactions 

Number  of  species 

Pressure,  N/m^ 

Heat  conduction,  Wlm^ 

Scalar  conservation  variable 
Vector  of  conserved  variables 
Position  vector,  m 
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R 

Re 

Rn 

R 


S 

S 

T 

r 

TV 

TVD 

r 


^char 

♦ 

M 


u 


u.v.w 

/vV 


V 


0) 


jf.y.z 

/ 

Y 

Z 


a 


O-y 

P 


Ratio  of  differences 
Gas  constant,  J/kg  K 
Reynolds  number 
Nose  radius,  m 

Universal  gas  constant,  Jlkmole  K 
Residual 

Generic  differences 

■j 

Surface  area,  m 
Matrix  of  eigenvectors 
Temperature,  K 
Integration  variable,  K 
Total  variation 
Total  variation  diminishing 
In  T 
Time,  s 

Characteristic  flow  time,  s 
Velocity  vector 
Diffusion  velocity  vector 
Cartesian  velocities,  mts 
Diffusion  velocities,  mis 
Cell  volume, 

Chemical  source  term,  kgln^lsec 
Chemical  source  term,  kmolelm^isec 
Generic  nonequilibrium  energy  source  term,  Jlm^/sec 
Cartesian  coordinates,  m 
y  ^  I  nondimensional  viscous  spacing 

Mass  fraction 
Partition  function 
Exponent  in  Park' s  TTy  model 
Viscous  relaxation  parameter 
dpIhE 
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T 

€ 

e 

K 

•S-N 

X 

A 

Ji 


V 


Vi 

n 

p 

0 


X 

¥ 

¥ 

<D 

J 

“l 
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Concentration,  kmoleln? 

Frozen  ratio  of  specific  heats 
Dissipation  parameter 
Energy  of  electronic  state,  J/kmole 
Computational  coordinates 
Characteristic  temperature,  K 
Thermal  conductivity,  Jims  K 
Dissipation  parameter 
Spectral  radius,  rr?ls 
Diagonal  matrix  of  eigenvalues 
Viscosity  of  mixture,  kglms 
Reduced  mass 
Pressure  smoothness  sensors 
Stoichiometric  coefficients 
Characteristic  vibrational  frequency 
Source  vector 
Density,  kgtrr^ 

Collision  cross  section, 

Directed  surface  area  at  cell  face 

Directed  surface  area  at  cell  center 

Relaxation  time,  s 

Shear  stress  tensor 

Flux  limiter 

Flux  limiter 

Mole  fraction 

Flux  limiter 

Dimensionless  electronic  energy 
Mass  source  of  species  s,  kg/m^/s 
Vibrational/electronic  source  term,  Jim  is 
Total  derivative 
Partial  derivative 


131 


AEDC-TR-94-18 


-» 

V 

E)el  operator 

V  (.)  " 

•»  J 

5 

Kronecker  delta 

6" 

Change  in  mole  number 

8(.)" 

J  ] 

-  w;.,« 

A(.)" 
r  V 

wr'  -(•>; 

A  (•)  '■ 

■*  0 

An 

linitial  viscous  spacing,  m 

52 

Implicit  time  change  in  vector  of  conserved  quantities 

AQ 

Explicit  time  change  in  vector  of  conserved  quantities 

A  f 

Time  step,  s 

Differences  in  computational  coordinates 

Subscripts: 

e 

Electronic 

f 

Frozen 

I 

Inviscid 

i,},  k 

Index  notation 

/ 

Interna]  mode 

L,R 

Left,  right 

n,m 

Summation/iteration  index 

0 

Stagnation/total  condition 

o,  e 

Odd/even  indices 

9 

Generic  temperature,  K 

r,  j 

Species/reaction  index 

r 

Rotational 

t 

Translational 

tr 

Translational-rotational 

V 

Viscous 

V 

Vibrational 

V 

Vibrational-electronic 
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V  -  c  Vibration-chemistry 

I 

v-T  Vibration-translational 

oo  Free-stream  condition 

Superscripts: 

atm  Standard  atmosphere 

/  Forward  rate 

H  Higher  order  approximation 

L  Lower  order  approximation 

L,  R  Left,  right 

n  Time  level 

o,  e  Odd/even  indices 

m  Iteration  level 

r  Reverse  rate 

s  Denotes  species  value 

V  Vibrational 

Equilibrium,  nozzle  irUet,  or  latest  value 
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